C SOURCES:
C    Simone Paziani, Saverio Moroni, Paola Gori-Giorgi, and Giovanni B. Bachelet.
C    Local-spin-densityfunctional for multideterminant density functional theory.
C    Physical Review B, 73(15), apr 2006.



C*****************************************************************************
      pure subroutine ESRC_VWN5_ERF(rho_c, mu, E)
C*****************************************************************************
C   Implemented by E.R. Kjellgren.
C
C   Subroutine generated using Sympy 1.3
C   Generated: March 25, 2019
C*****************************************************************************
      implicit none
      real*8, intent(in) :: rho_c, mu
      real*8, intent(out) :: E
      real*8 :: x0, x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12, x
     &13, x14, x15, x16, x17, x18, x19, x20, x21, x22, x23, x24, x25
      E = 0.0d0
      x0 = mu**2
      x1 = 2.14502939711103d0
      x2 = rho_c**0.666666666666667d0
      x3 = 1/(x1*x2)
      x4 = 1.41421356237310d0
      x5 = 0.682784063255296d0*rho_c**(-0.333333333333333d0)
      x6 = 4.60115111447049d0
      x7 = rho_c**(-1.33333333333333d0)/x6
      x8 = 3.14159265358979d0
      x9 = rho_c**1.0d0
      x10 = 1d0/x9
      x11 = x10/x8
      x12 = (-0.009578475d0*x11 + 0.0676322201645715d0*x3 + 0.0188171922
     &990732d0*x5 + 0.00126674656487366d0*x7 + 1.0d0)*exp(-0.68361076118
     &67116d0*x5)
      x13 = rho_c**(-1.66666666666667d0)
      x14 = 0.0837224526465878d0*x13
      x15 = x0*x5
      x16 = 0.826307487110758d0*rho_c**(-0.166666666666667d0)
      x17 = 5.57236303607474d0*mu*x16 + 1.0d0
      x18 = mu**3
      x19 = 0.90856029641607d0*x5
      x20 = 1d0/(3.5529372610885d0*x16 + x19 + 12.9352d0)
      x21 = 0.0310907d0*log(x19*x20) + 0.000969022771154437d0*log(x20*( 
     &0.95318429299693652d0*x16 + 0.10498d0)**2) + 0.038783294878113d0* 
     &atan(6.1519908197590798d0/(1.906368585993873d0*x16 + 3.72744000000
     &00001d0))
      x22 = 2.28942848510666d0*x6
      x23 = 1.63540853354893d0*x22*(-0.0259335011650457d0*x5 + 1.0d0)/( 
     &0.0524148278841779d0*x3 + 0.494402081358784d0*x5 + 1.0d0)
      x24 = 0.148394183600699d0
      x25 = x12 - 1.0d0
      E = rho_c*(x21 - (0.00316102962473761d0*mu**8*rho_c**(-2.666666666
     &66667d0)* x21 + mu**6*(0.0533250677421794d0*rho_c**(-2.0d0)*x21 - 
     &0.0167302167879722d0*x13*x24*x25) - 0.0223069557172962d0*mu**5* x1
     &2*x14*x4 - mu**4*(0.00139418473233101d0*rho_c**(-1.0d0)*x24*( 10.9
     &027235569928d0*x1*(0.558025705063192d0*x3 - 0.352521395009435d0*x5
     &)*exp(-0.49698248213959023d0*x5) - 1.63540853354893d0*x22 + x23) +
     & 0.13157433081915d0*x11*x25 - 1.55214407110551d0*x21*x7) - x18*x4*
     &(0.0315054072231411d0*x10*x12 + 0.00111534778586481d0*x14*(x2*x23 
     &+ 12.0d0*x8*x9*( 0.825481812223657d0*x3 - 4.49737346725955d0*x5)*e
     &xp( -0.28165369188898165d0*x5))) - 0.0621813817393098d0*log(( 1.91
     &40710242289742d0*rho_c**(-0.5d0)*x18 + 6.7683429898050367d0* x15 +
     & x17)/(3.1331792677937806d0*x15 + x17)))/( 0.508616435555896d0*x0*
     &x3 + 1.0d0)**4)
      end subroutine


C*****************************************************************************
      pure subroutine D1ESRC_VWN5_ERF(rho_c, mu, E, d1E)
C*****************************************************************************
C   Implemented by E.R. Kjellgren.
C
C   Subroutine generated using Sympy 1.3
C   Generated: March 25, 2019
C*****************************************************************************
      implicit none
      real*8, intent(in) :: rho_c, mu
      real*8, intent(out) :: E, d1E(9)
      real*8 :: x0, x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12, x
     &13, x14, x15, x16, x17, x18, x19, x20, x21, x22, x23, x24, x25, x2
     &6, x27, x28, x29, x30, x31, x32, x33, x34, x35, x36, x37, x38, x39
     &, x40, x41, x42, x43, x44, x45, x46, x47, x48, x49, x50, x51, x52,
     & x53, x54, x55, x56, x57, x58, x59, x60, x61, x62, x63, x64, x65, 
     &x66, x67, x68, x69, x70, x71, x72, x73, x74, x75, x76, x77, x78, x
     &79, x80, x81, x82, x83, x84, x85, x86, x87, x88, x89, x90, x91, x9
     &2, x93, x94, x95, x96, x97, x98, x99, x100, x101, x102, x103, x104
     &, x105, x106, x107, x108, x109, x110, x111, x112, x113, x114, x115
     &, x116, x117, x118
      E = 0.0d0
      d1E(:) = 0.0d0
      x0 = mu**2
      x1 = 2.14502939711103d0
      x2 = 1d0/x1
      x3 = rho_c**0.666666666666667d0
      x4 = 1d0/x3
      x5 = x2*x4
      x6 = 0.508616435555896d0*x0*x5 + 1.0d0
      x7 = x6**(-4)
      x8 = mu**4
      x9 = 1.46459188756152d0
      x10 = 1d0/x9
      x11 = rho_c**0.333333333333333d0
      x12 = 1d0/x11
      x13 = x10*x12
      x14 = exp(-0.6836107611867116d0*x13)
      x15 = rho_c**(-1.33333333333333d0)
      x16 = 4.60115111447049d0
      x17 = 1d0/x16
      x18 = x15*x17
      x19 = rho_c**1.0d0
      x20 = 1d0/x19
      x21 = 3.14159265358979d0
      x22 = 1d0/x21
      x23 = 0.009578475d0*x22
      x24 = 0.0188171922990732d0*x13 + 0.00126674656487366d0*x18 - x20*x
     &23 + 0.0676322201645715d0*x5 + 1.0d0
      x25 = x14*x24
      x26 = x25 - 1.0d0
      x27 = 0.13157433081915d0*x26
      x28 = x20*x22
      x29 = x27*x28
      x30 = 0.148394183600699d0
      x31 = 2.28942848510666d0
      x32 = x16*x31
      x33 = 1.63540853354893d0*x32
      x34 = exp(-0.49698248213959023d0*x13)
      x35 = x34*(-0.352521395009435d0*x13 + 0.558025705063192d0*x5)
      x36 = x1*x35
      x37 = 10.9027235569928d0*x36
      x38 = 0.494402081358784d0*x13 + 0.0524148278841779d0*x5 + 1.0d0
      x39 = 1d0/x38
      x40 = -0.0259335011650457d0*x13 + 1.0d0
      x41 = x39*x40
      x42 = 1.63540853354893d0*x32
      x43 = x41*x42
      x44 = 0.00139418473233101d0*rho_c**(-1.0d0)*x30*(-x33 + x37 + x43)
     &
      x45 = 0.826307487110758d0
      x46 = rho_c**(-0.166666666666667d0)*x45
      x47 = 1.90636858599387d0*x46 + 3.72744d0
      x48 = 0.90856029641607d0*x13
      x49 = 1d0/(3.5529372610885d0*x46 + x48 + 12.9352d0)
      x50 = 0.953184292996937d0*x46 + 0.10498d0
      x51 = 0.0310907d0*log(x48*x49) + 0.000969022771154437d0*log(x49*x5
     &0**2) + 0.038783294878113d0*atan(6.1519908197590798d0/x47)
      x52 = 1.55214407110551d0*x18
      x53 = x51*x52
      x54 = rho_c**(-1.66666666666667d0)
      x55 = x30*x54
      x56 = 0.101321183642338d0
      x57 = x51*x56
      x58 = rho_c**(-2.0d0)
      x59 = 0.526297323276602d0*x58
      x60 = 1.41421356237310d0
      x61 = mu**3
      x62 = 0.179587122125167d0
      x63 = 0.175432441092201d0*x20*x62
      x64 = 0.0837224526465878d0
      x65 = exp(-0.28165369188898165d0*x13)
      x66 = x65*(-4.49737346725955d0*x13 + 0.825481812223657d0*x5)
      x67 = 12.0d0*x21
      x68 = x66*x67
      x69 = x64*(x19*x68 + x3*x43)
      x70 = 0.00111534778586481d0*x54
      x71 = x0*x13
      x72 = 5.57236303607474d0*mu*x46 + 1.0d0
      x73 = 0.564189583547756d0
      x74 = 3.39260255780131d0*rho_c**(-0.5d0)*x61*x73 + 6.7683429898050
     &4d0*x71 + x72
      x75 = x74/(3.13317926779378d0*x71 + x72)
      x76 = 0.0621813817393098d0
      x77 = 0.0472353356922751d0
      x78 = x51*x77
      x79 = rho_c**(-2.66666666666667d0)
      x80 = 0.0669208671518886d0*x79
      x81 = mu**5
      x82 = x25*x60
      x83 = x64*x82
      x84 = 0.0223069557172962d0*x54
      x85 = mu**8*x78*x80 + mu**6*(-0.0167302167879722d0*x26*x55 + x57*x
     &59) - x60* x61*(x25*x63 + x69*x70) - x76*log(x75) - x81*x83*x84
      x86 = x51 - x7*(-x8*(x29 + x44 - x53) + x85)
      x87 = x47**(-2)
      x88 = rho_c**(-1.16666666666667d0)*x45
      x89 = 0.0758081683534927d0*x87*x88/(37.8469910464d0*x87 + 1.0d0)
      x90 = 0.30285343213869d0*x15
      x91 = x49*(x10*x90 + 0.592156210181417d0*x88)
      x92 = 0.90856029641607d0*x12*x91
      x93 = 0.0342197431724027d0*x11
      x94 = x93*(x90 - x92)
      x95 = 1d0/x50
      x96 = 0.000307885761673592d0*x88
      x97 = 0.000969022771154437d0*x50*x91
      x98 = x95*(x96 - x97)
      x99 = x2*x54
      x100 = rho_c**(-2.33333333333333d0)
      x101 = x100*x17
      x102 = 0.00168899541983155d0*x101
      x103 = x23*x58
      x104 = 0.045088146776381d0*x99
      x105 = x10*x15
      x106 = 0.00627239743302441d0*x105
      x107 = rho_c**(-3.0d0)
      x108 = mu**7
      x109 = 0.928727172679123d0*x88
      x110 = mu*x105
      x111 = x22*x58
      x112 = x14*(-0.000844497709915773d0*x101 + 0.113935126864452d0*x10
     &5*x24 - 0.0031361987165122d0*x105 + 0.0047892375d0*x111 - 0.022544
     &0733881905d0*x99)
      x113 = x89 + x93*(-x90 + x92) + x95*(-x96 + x97)
      x114 = 0.0141372897033723d0*x21*x31*x39*x4
      x115 = rho_c**(-0.333333333333333d0)
      x116 = x115*x32
      x117 = 1.09027235569928d0*x116*x41
      x118 = x3*x40*x42*(0.164800693786261d0*x105 + 0.034943218589452d0*
     &x99)/x38**2
      E = rho_c*x86
      d1E(1) = -rho_c*(mu*x7*(-0.17845564573837d0*rho_c**(-3.66666666666
     &667d0)*x108*x78 - x0*x60*(0.00490180588786583d0*x100*x25 + x14*x63
     &*(-x102 + x103 - x104 - x106) - 0.175432441092201d0*x25*x58*x62 + 
     &x64*x70*( 2.41662179562687d0*rho_c**(-0.333333333333333d0)*x66 + x
     &114 + x117 + x118 + x19*x65*x67*(1.49912448908652d0*x105 - 0.55032
     &1208149104d0*x99) + x68) - 0.00185891297644135d0*x69*x79) - 0.0002
     &90571663240495d0*x107*x8*x82 + x108*x77*x80*(x89 - x94 - x98) + x1
     &4*x60*x64*x8*x84*(x102 - x103 + x104 + x106) - x61*( 2.06952542814
     &068d0*x101*x51 - x111*x27 + 0.263148661638301d0*x112 *x28 - x113*x
     &52 + 0.00232364122055169d0*x30*x58*(x33 - x37 - x43 ) - 0.00139418
     &473233101d0*x55*(-10.9027235569928d0*x1*x3*x34*( 0.117507131669812
     &d0*x105 - 0.372017136708795d0*x99) - x114 - 7.26848237132856d0*x11
     &5*x36 + 1.09027235569928d0*x116 - x117 - x118 - 1.80615420514536d0
     &*x35*x4*x9)) + 0.037178259528827d0*x79* x8*x83 - x81*(1.0525946465
     &532d0*x107*x57 + 0.0334604335759443d0* x112*x55 - x113*x56*x59 - 0
     &.0278836946466203d0*x26*x30*x79) - x76 *(-1.69630127890066d0*rho_c
     &**(-1.5d0)*x0*x73 - x109 - 2.25611432993501d0*x110 + x75*(x109 + 1
     &.04439308926459d0*x110))/ x74) + 1.35631049481572d0*x0*x99*(x8*(-x
     &29 - x44 + x53) + x85)/x6 **5 - x89 + x94 + x98) + x86
      end subroutine


C*****************************************************************************
      pure subroutine D2ESRC_VWN5_ERF(rho_c, mu, E, d1E, d2E)
C*****************************************************************************
C   Implemented by E.R. Kjellgren.
C
C   Subroutine generated using Sympy 1.3
C   Generated: March 25, 2019
C*****************************************************************************
      implicit none
      real*8, intent(in) :: rho_c, mu
      real*8, intent(out) :: E, d1E(9), d2E(45)
      real*8 :: x0, x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12, x
     &13, x14, x15, x16, x17, x18, x19, x20, x21, x22, x23, x24, x25, x2
     &6, x27, x28, x29, x30, x31, x32, x33, x34, x35, x36, x37, x38, x39
     &, x40, x41, x42, x43, x44, x45, x46, x47, x48, x49, x50, x51, x52,
     & x53, x54, x55, x56, x57, x58, x59, x60, x61, x62, x63, x64, x65, 
     &x66, x67, x68, x69, x70, x71, x72, x73, x74, x75, x76, x77, x78, x
     &79, x80, x81, x82, x83, x84, x85, x86, x87, x88, x89, x90, x91, x9
     &2, x93, x94, x95, x96, x97, x98, x99, x100, x101, x102, x103, x104
     &, x105, x106, x107, x108, x109, x110, x111, x112, x113, x114, x115
     &, x116, x117, x118, x119, x120, x121, x122, x123, x124, x125, x126
     &, x127, x128, x129, x130, x131, x132, x133, x134, x135, x136, x137
     &, x138, x139, x140, x141, x142, x143, x144, x145, x146, x147, x148
     &, x149, x150, x151, x152, x153, x154, x155, x156, x157, x158, x159
     &, x160, x161, x162, x163, x164, x165, x166, x167, x168, x169, x170
     &, x171, x172, x173, x174, x175, x176, x177, x178, x179, x180, x181
     &, x182, x183, x184, x185, x186, x187, x188, x189, x190, x191, x192
     &, x193, x194, x195, x196, x197, x198, x199, x200, x201, x202, x203
     &, x204, x205, x206, x207, x208, x209, x210, x211, x212, x213, x214
     &, x215, x216, x217, x218, x219, x220, x221, x222, x223, x224, x225
     &, x226, x227, x228, x229, x230, x231, x232, x233, x234, x235, x236
     &, x237, x238, x239, x240, x241, x242, x243, x244, x245, x246, x247
     &, x248
      E = 0.0d0
      d1E(:) = 0.0d0
      d2E(:) = 0.0d0
      x0 = mu**2
      x1 = 2.14502939711103d0
      x2 = 1d0/x1
      x3 = rho_c**0.666666666666667d0
      x4 = 1d0/x3
      x5 = x2*x4
      x6 = 0.508616435555896d0*x0*x5 + 1.0d0
      x7 = x6**(-4)
      x8 = mu**4
      x9 = 1.46459188756152d0
      x10 = 1d0/x9
      x11 = rho_c**0.333333333333333d0
      x12 = 1d0/x11
      x13 = x10*x12
      x14 = exp(-0.6836107611867116d0*x13)
      x15 = rho_c**(-1.33333333333333d0)
      x16 = 4.60115111447049d0
      x17 = 1d0/x16
      x18 = x15*x17
      x19 = rho_c**1.0d0
      x20 = 1d0/x19
      x21 = 3.14159265358979d0
      x22 = 1d0/x21
      x23 = 0.009578475d0*x22
      x24 = 0.0188171922990732d0*x13 + 0.00126674656487366d0*x18 - x20*x
     &23 + 0.0676322201645715d0*x5 + 1.0d0
      x25 = x14*x24
      x26 = x25 - 1.0d0
      x27 = 0.13157433081915d0*x26
      x28 = x20*x22
      x29 = x27*x28
      x30 = 2.28942848510666d0
      x31 = x16*x30
      x32 = 1.63540853354893d0*x31
      x33 = exp(-0.49698248213959023d0*x13)
      x34 = x33*(-0.352521395009435d0*x13 + 0.558025705063192d0*x5)
      x35 = x1*x34
      x36 = 10.9027235569928d0*x35
      x37 = 0.494402081358784d0*x13 + 0.0524148278841779d0*x5 + 1.0d0
      x38 = 1d0/x37
      x39 = -0.0259335011650457d0*x13 + 1.0d0
      x40 = x38*x39
      x41 = 1.63540853354893d0*x31
      x42 = x40*x41
      x43 = 0.148394183600699d0
      x44 = 0.00139418473233101d0*rho_c**(-1.0d0)*x43
      x45 = x44*(-x32 + x36 + x42)
      x46 = 0.826307487110758d0
      x47 = rho_c**(-0.166666666666667d0)*x46
      x48 = 1.90636858599387d0*x47 + 3.72744d0
      x49 = 0.90856029641607d0*x13
      x50 = 3.5529372610885d0*x47 + x49 + 12.9352d0
      x51 = 1d0/x50
      x52 = 0.953184292996937d0*x47 + 0.10498d0
      x53 = x52**2
      x54 = x51*x53
      x55 = 0.000969022771154437d0*log(x54) + 0.0310907d0*log(x49*x51) +
     & 0.038783294878113d0*atan(6.1519908197590798d0/x48)
      x56 = 1.55214407110551d0*x18
      x57 = x55*x56
      x58 = 0.693147180559945d0
      x59 = -x58 + 1.0d0
      x60 = x0*x13
      x61 = 5.57236303607474d0*mu*x47 + 1.0d0
      x62 = 3.13317926779378d0*x60 + x61
      x63 = 1d0/x62
      x64 = 0.564189583547756d0
      x65 = mu**3
      x66 = 3.39260255780131d0*rho_c**(-0.5d0)*x64*x65 + 6.7683429898050
     &4d0*x60 + x61
      x67 = x63*x66
      x68 = 0.202642367284676d0
      x69 = x68*log(x67)
      x70 = rho_c**(-1.66666666666667d0)
      x71 = x43*x70
      x72 = 0.101321183642338d0
      x73 = x55*x72
      x74 = rho_c**(-2.0d0)
      x75 = 0.526297323276602d0*x74
      x76 = 1.41421356237310d0
      x77 = 0.179587122125167d0
      x78 = x25*x77
      x79 = 0.175432441092201d0*x20
      x80 = 0.0837224526465878d0
      x81 = -4.49737346725955d0*x13 + 0.825481812223657d0*x5
      x82 = exp(-0.28165369188898165d0*x13)
      x83 = x21*x82
      x84 = 12.0d0*x83
      x85 = x81*x84
      x86 = x80*(x19*x85 + x3*x42)
      x87 = 0.00111534778586481d0*x70
      x88 = 0.0472353356922751d0
      x89 = x55*x88
      x90 = rho_c**(-2.66666666666667d0)
      x91 = 0.0669208671518886d0*x90
      x92 = mu**5
      x93 = x25*x80
      x94 = x76*x93
      x95 = 0.0223069557172962d0*x70
      x96 = mu**8*x89*x91 + mu**6*(-0.0167302167879722d0*x26*x71 + x73*x
     &75) - x65* x76*(x78*x79 + x86*x87) - x92*x94*x95
      x97 = -x59*x69 + x96
      x98 = x55 - x7*(-x8*(x29 + x45 - x57) + x97)
      x99 = rho_c**(-1.16666666666667d0)*x46
      x100 = x48**(-2)
      x101 = 37.8469910464d0*x100 + 1.0d0
      x102 = 1d0/x101
      x103 = x100*x102
      x104 = x103*x99
      x105 = 0.0758081683534927d0*x104
      x106 = 0.30285343213869d0*x15
      x107 = x10*x106 + 0.592156210181417d0*x99
      x108 = x107*x51
      x109 = 0.90856029641607d0*x12
      x110 = x108*x109
      x111 = x106 - x110
      x112 = 0.0342197431724027d0*x11
      x113 = x111*x112
      x114 = 1d0/x52
      x115 = 0.000307885761673592d0*x99
      x116 = x108*x52
      x117 = 0.000969022771154437d0*x116
      x118 = x114*(x115 - x117)
      x119 = -x29 + x57
      x120 = x6**(-5)
      x121 = x2*x70
      x122 = x120*x121
      x123 = x0*x122*(x8*(x119 - x45) + x97)
      x124 = x22*x74
      x125 = x124*x27
      x126 = x32 - x36 - x42
      x127 = x126*x43
      x128 = 0.00232364122055169d0*x127*x74
      x129 = x10*x15
      x130 = rho_c**(-2.33333333333333d0)
      x131 = x130*x17
      x132 = -0.0225440733881905d0*x121 + 0.0047892375d0*x124 - 0.003136
     &1987165122d0* x129 - 0.000844497709915773d0*x131
      x133 = x14*(0.113935126864452d0*x129*x24 + x132)
      x134 = 0.263148661638301d0*x28
      x135 = x133*x134
      x136 = 2.06952542814068d0*x131*x55
      x137 = -x106 + x110
      x138 = x112*x137
      x139 = x105 + x114*(-x115 + x117) + x138
      x140 = x139*x56
      x141 = rho_c**(-0.333333333333333d0)
      x142 = x141*x31
      x143 = x4*x9
      x144 = x33*(-0.372017136708795d0*x121 + 0.117507131669812d0*x129)
      x145 = x1*x144
      x146 = 10.9027235569928d0*x3
      x147 = x21*x30*x4
      x148 = 0.0141372897033723d0*x147*x38
      x149 = 1.09027235569928d0*x142*x40
      x150 = x37**(-2)
      x151 = 0.034943218589452d0*x121 + 0.164800693786261d0*x129
      x152 = x150*x151
      x153 = x152*x39
      x154 = x3*x41
      x155 = x153*x154
      x156 = -7.26848237132856d0*x141*x35 + 1.09027235569928d0*x142 - 1.
     &80615420514536d0*x143*x34 - x145*x146 - x148 - x149 - x155
      x157 = 0.00139418473233101d0*x71
      x158 = x156*x157
      x159 = x43*x90
      x160 = 0.0278836946466203d0*x159*x26
      x161 = 0.0334604335759443d0*x71
      x162 = x133*x161
      x163 = rho_c**(-3.0d0)
      x164 = 1.0525946465532d0*x163*x73
      x165 = x72*x75
      x166 = x139*x165
      x167 = 0.00168899541983155d0*x131
      x168 = x23*x74
      x169 = 0.045088146776381d0*x121
      x170 = 0.00627239743302441d0*x129
      x171 = x14*(-x167 + x168 - x169 - x170)
      x172 = x171*x77
      x173 = x172*x79
      x174 = 0.122619224952946d0
      x175 = x174*x25
      x176 = 0.0399758348639607d0*x130*x175
      x177 = 0.175432441092201d0*x74*x78
      x178 = 0.00185891297644135d0*x86*x90
      x179 = x81*x82
      x180 = 2.14502939711103d0
      x181 = rho_c**(-0.333333333333333d0)*x180
      x182 = -0.550321208149104d0*x121 + 1.49912448908652d0*x129
      x183 = x19*x84
      x184 = x80*(x148 + x149 + x155 + 1.12661476755593d0*x179*x181 + x1
     &82*x183 + x85 )
      x185 = x184*x87
      x186 = x0*x76
      x187 = 0.928727172679123d0*x99
      x188 = mu*x129
      x189 = 2.25611432993501d0*x188
      x190 = x0*x64
      x191 = 1.69630127890066d0*rho_c**(-1.5d0)*x190
      x192 = x187 + 1.04439308926459d0*x188
      x193 = -x187 - x189 - x191 + x192*x67
      x194 = x68/x66
      x195 = x193*x194
      x196 = mu**7
      x197 = x196*x88
      x198 = x197*x91
      x199 = x76*x8
      x200 = x199*x80
      x201 = x200*x95
      x202 = rho_c**(-3.66666666666667d0)
      x203 = x196*x89
      x204 = x8*x94
      x205 = 0.0571643564037363d0
      x206 = x199*x25
      x207 = x205*x206
      x208 = -0.00508309165921971d0*x163*x207 - 0.17845564573837d0*x202*
     &x203 + 0.037178259528827d0*x204*x90
      x209 = mu*x7
      x210 = x209*(x14*x201*(x167 - x168 + x169 + x170) - x186*(x173 + x
     &176 - x177 - x178 + x185) - x195*x59 + x198*(x105 - x113 - x118) +
     & x208 - x65* (-x125 + x128 + x135 + x136 - x140 - x158) - x92*(-x1
     &60 + x162 + x164 - x166))
      x211 = 0.000615771523347183d0*x99
      x212 = rho_c**(-2.33333333333333d0)*x10
      x213 = 1.82319440383792d0*x212/(x101**2*x48**5)
      x214 = 0.0481727702369445d0*x102*x212/x48**3
      x215 = rho_c**(-2.16666666666667d0)*x46
      x216 = 0.0884428630790749d0*x103*x215
      x217 = 0.0114065810574676d0*rho_c**(-0.666666666666667d0)*x137
      x218 = 1d0/x53
      x219 = x116 - 0.317728097665645d0*x99
      x220 = x115*x218*x219
      x221 = x108*x138
      x222 = 0.000969022771154437d0*x108*x114*x219
      x223 = 0.40380457618492d0*x130
      x224 = x10*x223 + 0.69084891187832d0*x215
      x225 = x107*(0.60570686427738d0*x129 + 1.18431242036283d0*x99)/x50
     &**2
      x226 = x112*(-0.60570686427738d0*x108*x15 - x109*x224*x51 + x109*x
     &225 + x223)
      x227 = x218*(-x116*x211 + 4.89119786774443d-5*x212 + 0.00035920005
     &5285857d0* x215*x52 - 0.000969022771154437d0*x224*x54 + 0.00096902
     &2771154437d0*x225*x53)
      x228 = x58 - 1.0d0
      x229 = x228*x69 + x8*(x119 + x126*x44) + x96
      x230 = rho_c**(-3.33333333333333d0)
      x231 = x17*x230
      x232 = x2*x90
      x233 = x195*x228
      x234 = x163*x22
      x235 = x10*x130
      x236 = x14*(0.00394098931294027d0*x231 + 0.0751469112939684d0*x232
     & - 0.01915695d0*x234 + 0.00836319657736588d0*x235)
      x237 = rho_c**(-4.0d0)
      x238 = x187 + x189 + x191
      x239 = mu*x192
      x240 = x239*x63
      x241 = 1.08351503479231d0*x215
      x242 = mu*x235
      x243 = -x213 + x214 - x216 + x217 + x220 - x221 - x222 + x226 + x2
     &27
      x244 = rho_c**(-1.33333333333333d0)
      x245 = x244*x31
      x246 = x154*x39
      x247 = -2.18054471139857d0*x142*x153 - 0.0282745794067445d0*x147*x
     &152 + x150* x246*(0.0582386976490866d0*x232 + 0.219734258381682d0*
     &x235) - x151*x246*(0.0698864371789039d0*x121 + 0.329601387572523d0
     &*x129)/ x37**3 + 0.363424118566428d0*x245*x40
      x248 = x14*(0.455740507457808d0*x129*x132 - x163*x23 + 0.001970494
     &65647014d0* x231 + 0.0259624262672375d0*x232*x24 + 0.0375734556469
     &842d0*x232 - 0.151913502485936d0*x235*x24 + 0.00418159828868294d0*
     &x235)
      E = rho_c*x98
      d1E(1) = -rho_c*(-x105 + x113 + x118 + 1.35631049481572d0*x123 + x
     &210) + x98
      d2E(1) = -rho_c*(-2.26051749135954d0*x0*x120*x229*x232 + 2.7126209
     &8963145d0*x122* x65*(x139*x198 - x171*x201 + x186*(-x173 - x176 + 
     &x177 + x178 - x185) + x208 + x233 + x65*(x125 - x128 - x135 - x136
     & + x140 + x158) + x92*(x160 - x162 - x164 + x166)) + x209*(mu*x193
     &*x228* x238*x68/x66**2 + 0.654337367707355d0*rho_c**(-4.6666666666
     &6667d0 )*x203 - 4.52089344419912d-5*rho_c**(-4.33333333333333d0)*x
     &206 - 0.356911291476739d0*x139*x197*x202 - 0.0101661833184394d0*x1
     &63* x171*x199*x205 + 0.074356519057654d0*x171*x200*x90 + x186*( -0
     &.00910930363347549d0*rho_c**(-3.66666666666667d0)*x93 - 0.07995166
     &97279215d0*x130*x171*x174 - 0.350864882184401d0*x163* x78 + 0.3508
     &64882184401d0*x172*x74 + 0.133252782879869d0*x175* x230 + 0.003717
     &8259528827d0*x184*x90 - 0.00495710127051027d0*x202 *x86 - x236*x77
     &*x79 + x80*x87*(-0.751076511703951d0*x15*x179*x180 - 0.15491242678
     &0983d0*x179*x70 - 2.25322953511185d0*x181*x182*x82 - 24.0d0*x182*x
     &83 - x183*(0.917202013581841d0*x232 - 1.99883265211535d0*x235) + x
     &247)) + x194*x228*(2.54445191835098d0 *rho_c**(-2.5d0)*x190 - 2.0d
     &0*x238*x240 + x239*x66*( 2.08878617852919d0*x188 + 1.8574543453582
     &5d0*x99)/x62**2 + x241 + 3.00815243991335d0*x242 - x67*(x241 + 1.3
     &9252411901946d0*x242)) + x198*x243 - x201*x236 - 0.099142025410205
     &4d0*x202*x204 + 0.023721094409692d0*x207*x237 - x233*x240 + x65*( 
     &0.526297323276602d0*x124*x133 + 0.00619637658813783d0*x127*x163 - 
     &4.13905085628136d0*x131*x139 - x134*x248 - 0.00464728244110338d0* 
     &x156*x159 + x157*(-x1*x146*x33*(0.620028561181324d0*x232 - 0.15667
     &6175559749d0*x235) - 14.5369647426571d0*x141*x145 - 3.612308410290
     &72d0*x143*x144 + 2.42282745710952d0*x244*x35 - 0.363424118566428d0
     &*x245 + x247 - 0.299209d0*x34*x74) + 4.82889266566159d0*x231*x55 -
     & 0.263148661638301d0*x234*x26 + x243 *x56) + x92*(0.11153477858648
     &1d0*x133*x159 - 2.10518929310641d0* x139*x163*x72 - x161*x248 + x1
     &65*x243 - 0.074356519057654d0*x202* x26*x43 + 3.15778393965961d0*x
     &237*x73)) + x213 - x214 + x216 - x217 - x220 + x221 + x222 - x226 
     &- x227 + 2.29947269793409d0*x229 *x231*x8/x6**6) + 0.1516163367069
     &85d0*x104 - 0.0684394863448054d0 *x11*x111 - x114*(-0.001938045542
     &30887d0*x116 + x211) - 2.71262098963145d0*x123 - 2.0d0*x210
      end subroutine


C*****************************************************************************
      pure subroutine ESRC_SPIN_VWN5_ERF(rho_c, rho_s, mu, E)
C*****************************************************************************
C   Implemented by E.R. Kjellgren.
C
C   Subroutine generated using Sympy 1.3
C   Generated: March 25, 2019
C*****************************************************************************
      implicit none
      real*8, intent(in) :: rho_c, rho_s, mu
      real*8, intent(out) :: E
      real*8 :: x0, x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12, x
     &13, x14, x15, x16, x17, x18, x19, x20, x21, x22, x23, x24, x25, x2
     &6, x27, x28, x29, x30, x31, x32, x33, x34, x35, x36, x37, x38, x39
     &, x40, x41, x42, x43, x44, x45, x46, x47, x48, x49, x50, x51, x52,
     & x53, x54
      E = 0.0d0
      x0 = 0.826307487110758d0*rho_c**(-0.166666666666667d0)
      x1 = 1.90636858599387d0*x0
      x2 = 0.038783294878113d0*atan(6.1519908197590798d0/(x1 + 3.7274400
     &000000001d0 ))
      x3 = 0.682784063255296d0*rho_c**(-0.333333333333333d0)
      x4 = 0.90856029641607d0*x3
      x5 = 1d0/(3.5529372610885d0*x0 + x4 + 12.9352d0)
      x6 = 0.0310907d0*log(x4*x5)
      x7 = 0.953184292996937d0*x0
      x8 = 0.000969022771154437d0*log(x5*(x7 + 0.10498d0)**2)
      x9 = 0.101321183642338d0
      x10 = rho_s**4/rho_c**4
      x11 = rho_s/rho_c
      x12 = x11 + 1.0d0
      x13 = -x11 + 1.0d0
      x14 = 1.92366105093154d0*x12**1.33333333333333d0 + 1.9236610509315
     &4d0*x13** 1.33333333333333d0 - 3.84732210186307d0
      x15 = 1d0/(1.07811815828004d0*x0 + x4 + 13.0045d0)
      x16 = 0.0974703937105774d0*x14*x9*(-x10 + 1.0d0)*(log(x15*x4) + 0.
     &000414033794282063d0*log(x15*(x7 + 0.0047584000000000003d0)**2 ) +
     & 0.317708004743941d0*atan(7.1231089178181177d0/(x1 + 1.13107d0 )))
     &
      x17 = 1d0/(6.72988144596143d0*x0 + x4 + 18.0578d0)
      x18 = x2 + x6 + x8
      x19 = x10*x14*(x18 - 0.01554535d0*log(x17*x4) - 0.0022478670955426
     &1d0*log(x17* (x7 + 0.32500000000000001d0)**2) - 0.0524913931697809
     &d0*atan( 4.7309269095601136d0/(x1 + 7.0604199999999997d0)))
      x20 = mu**2
      x21 = 2.14502939711103d0
      x22 = rho_c**0.666666666666667d0
      x23 = 1/(x21*x22)
      x24 = 1.41421356237310d0
      x25 = rho_c**(-1.66666666666667d0)
      x26 = 0.0837224526465878d0*x25
      x27 = 1.0d0 - rho_s**2/rho_c**2
      x28 = 4.60115111447049d0
      x29 = rho_c**(-1.33333333333333d0)/x28
      x30 = 3.14159265358979d0
      x31 = rho_c**1.0d0
      x32 = 1d0/x31
      x33 = x32/x30
      x34 = (0.0676322201645715d0*x23 + 0.00126674656487366d0*x29 + 0.01
     &88171922990732d0*x3 - 0.009578475d0*x33 + 1.0d0)*exp( -0.683610761
     &1867116d0*x3)
      x35 = x27*x34
      x36 = x12**0.666666666666667d0 + x13**0.666666666666667d0
      x37 = x36**3
      x38 = x20*x3/x36**2
      x39 = 11.1447260721495d0*mu*x0/x36 + 1.0d0
      x40 = mu**3
      x41 = 1d0/x12
      x42 = x41**0.666666666666667d0
      x43 = 0.0524148278841779d0*x23
      x44 = x3*x41**0.333333333333333d0
      x45 = 2.28942848510666d0*x28
      x46 = 0.817704266774463d0*x45
      x47 = x12**2*x46*(-0.0259335011650457d0*x44 + 1.0d0)/(x42*(x42*x43
     & + 0.494402081358784d0*x44 + 1.0d0))
      x48 = 1d0/x13
      x49 = x48**0.666666666666667d0
      x50 = x3*x48**0.333333333333333d0
      x51 = x13**2*x46*(-0.0259335011650457d0*x50 + 1.0d0)/(x49*(x43*x49
     & + 0.494402081358784d0*x50 + 1.0d0))
      x52 = -x16 + x18 - x19
      x53 = 0.148394183600699d0
      x54 = x27*(x34 - 1.0d0)
      E = -rho_c*(x16 + x19 - x2 - x6 - x8 + 0.25d0*(0.0126441184989504d
     &0*mu**8* rho_c**(-2.66666666666667d0)*x52 + 4.0d0*mu**6*( 0.053325
     &0677421794d0*rho_c**(-2.0d0)*x52 - 0.0167302167879722d0* x25*x53*x
     &54) - 0.0892278228691848d0*mu**5*x24*x26*x35 + 4.0d0*mu **4*(-0.00
     &139418473233101d0*rho_c**(-1.0d0)*x53*( 10.9027235569928d0*x21*x27
     &*(0.558025705063192d0*x23 - 0.352521395009435d0*x3)*exp(-0.4969824
     &8213959023d0*x3) - 0.817704266774463d0*x45*(x12**2.66666666666667d
     &0 + x13** 2.66666666666667d0) + x47 + x51) + 1.55214407110551d0*x2
     &9*x52 - 0.13157433081915d0*x33*x54) - 4.0d0*x24*x40*( 0.0011153477
     &8586481d0*x26*(x22*x47 + x22*x51 + 12.0d0*x27*x30*x31 *(0.82548181
     &2223657d0*x23 - 4.49737346725955d0*x3)*exp( -0.28165369188898165d0
     &*x3)) + 0.0315054072231411d0*x32*x35) - 0.306852819440055d0*x37*x9
     &*log((15.312568193831794d0*rho_c**( -0.5d0)*x40/x37 + 27.073371959
     &220147d0*x38 + x39)/( 12.532717071175123d0*x38 + x39)))/(0.5086164
     &35555896d0*x20*x23 + 1.0d0)**4)
      end subroutine


C*****************************************************************************
      pure subroutine D1ESRC_SPIN_VWN5_ERF(rho_c, rho_s, mu, E, d1E)
C*****************************************************************************
C   Implemented by E.R. Kjellgren.
C
C   Subroutine generated using Sympy 1.3
C   Generated: March 25, 2019
C*****************************************************************************
      implicit none
      real*8, intent(in) :: rho_c, rho_s, mu
      real*8, intent(out) :: E, d1E(9)
      real*8 :: x0, x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12, x
     &13, x14, x15, x16, x17, x18, x19, x20, x21, x22, x23, x24, x25, x2
     &6, x27, x28, x29, x30, x31, x32, x33, x34, x35, x36, x37, x38, x39
     &, x40, x41, x42, x43, x44, x45, x46, x47, x48, x49, x50, x51, x52,
     & x53, x54, x55, x56, x57, x58, x59, x60, x61, x62, x63, x64, x65, 
     &x66, x67, x68, x69, x70, x71, x72, x73, x74, x75, x76, x77, x78, x
     &79, x80, x81, x82, x83, x84, x85, x86, x87, x88, x89, x90, x91, x9
     &2, x93, x94, x95, x96, x97, x98, x99, x100, x101, x102, x103, x104
     &, x105, x106, x107, x108, x109, x110, x111, x112, x113, x114, x115
     &, x116, x117, x118, x119, x120, x121, x122, x123, x124, x125, x126
     &, x127, x128, x129, x130, x131, x132, x133, x134, x135, x136, x137
     &, x138, x139, x140, x141, x142, x143, x144, x145, x146, x147, x148
     &, x149, x150, x151, x152, x153, x154, x155, x156, x157, x158, x159
     &, x160, x161, x162, x163, x164, x165, x166, x167, x168, x169, x170
     &, x171, x172, x173, x174, x175, x176, x177, x178, x179, x180, x181
     &, x182, x183, x184, x185, x186, x187, x188, x189, x190, x191, x192
     &, x193, x194, x195, x196, x197, x198, x199, x200, x201, x202, x203
     &, x204, x205, x206, x207, x208, x209, x210, x211, x212, x213, x214
     &, x215, x216, x217, x218, x219, x220, x221, x222, x223, x224, x225
     &, x226, x227, x228, x229, x230, x231, x232, x233, x234, x235, x236
     &, x237, x238, x239, x240, x241, x242, x243, x244, x245, x246, x247
     &, x248, x249, x250, x251, x252, x253, x254, x255, x256, x257, x258
     &, x259, x260, x261, x262, x263, x264, x265, x266, x267, x268, x269
     &, x270, x271, x272, x273, x274, x275, x276, x277, x278, x279, x280
     &
      E = 0.0d0
      d1E(:) = 0.0d0
      x0 = rho_c**(-2)
      x1 = rho_s**2
      x2 = x0*x1
      x3 = -x2 + 1.0d0
      x4 = 1.46459188756152d0
      x5 = 1d0/x4
      x6 = rho_c**0.333333333333333d0
      x7 = 1d0/x6
      x8 = x5*x7
      x9 = exp(-0.6836107611867116d0*x8)
      x10 = rho_c**(-1.33333333333333d0)
      x11 = 4.60115111447049d0
      x12 = 1d0/x11
      x13 = x10*x12
      x14 = rho_c**1.0d0
      x15 = 1d0/x14
      x16 = 3.14159265358979d0
      x17 = 1d0/x16
      x18 = 0.009578475d0*x17
      x19 = 2.14502939711103d0
      x20 = 1d0/x19
      x21 = rho_c**0.666666666666667d0
      x22 = 1d0/x21
      x23 = x20*x22
      x24 = 0.00126674656487366d0*x13 - x15*x18 + 0.0676322201645715d0*x
     &23 + 0.0188171922990732d0*x8 + 1.0d0
      x25 = x24*x9
      x26 = 0.0837224526465878d0
      x27 = 1.41421356237310d0
      x28 = mu**5*x27
      x29 = x26*x28
      x30 = x25*x29
      x31 = rho_c**(-1.66666666666667d0)
      x32 = 0.0892278228691848d0*x31
      x33 = 0.101321183642338d0
      x34 = 0.693147180559945d0
      x35 = x33*(-x34 + 1.0d0)
      x36 = 1d0/rho_c
      x37 = rho_s*x36
      x38 = x37 + 1.0d0
      x39 = -x37 + 1.0d0
      x40 = x38**0.666666666666667d0 + x39**0.666666666666667d0
      x41 = x40**3
      x42 = mu**2
      x43 = x40**2
      x44 = 1d0/x43
      x45 = x42*x44
      x46 = x45*x8
      x47 = 0.826307487110758d0
      x48 = rho_c**(-0.166666666666667d0)*x47
      x49 = 1d0/x40
      x50 = mu*x49
      x51 = 11.1447260721495d0*x48*x50 + 1.0d0
      x52 = mu**3
      x53 = 0.564189583547756d0
      x54 = x53/x41
      x55 = 27.1408204624105d0*rho_c**(-0.5d0)*x52*x54 + 27.073371959220
     &1d0*x46 + x51
      x56 = x55/(12.5327170711751d0*x46 + x51)
      x57 = log(x56)
      x58 = x41*x57
      x59 = 0.179587122125167d0
      x60 = x25*x59
      x61 = 0.175432441092201d0*x15
      x62 = x60*x61
      x63 = 0.825481812223657d0*x23 - 4.49737346725955d0*x8
      x64 = exp(-0.28165369188898165d0*x8)
      x65 = x16*x63*x64
      x66 = 12.0d0*x3
      x67 = x65*x66
      x68 = 2.28942848510666d0
      x69 = x11*x68
      x70 = 0.817704266774463d0*x69
      x71 = x38**2
      x72 = 1d0/x38
      x73 = x72**0.666666666666667d0
      x74 = 0.0524148278841779d0*x23
      x75 = x72**0.333333333333333d0
      x76 = x75*x8
      x77 = x73*x74 + 0.494402081358784d0*x76 + 1.0d0
      x78 = 1d0/x77
      x79 = x71*x78
      x80 = 1d0/x73
      x81 = 0.0259335011650457d0*x76
      x82 = x80*(-x81 + 1.0d0)
      x83 = x79*x82
      x84 = x70*x83
      x85 = 1d0/x39
      x86 = x85**0.666666666666667d0
      x87 = 1d0/x86
      x88 = x70*x87
      x89 = x85**0.333333333333333d0
      x90 = x8*x89
      x91 = 0.0259335011650457d0*x90
      x92 = -x91 + 1.0d0
      x93 = x39**2
      x94 = x74*x86 + 0.494402081358784d0*x90 + 1.0d0
      x95 = 1d0/x94
      x96 = x93*x95
      x97 = x92*x96
      x98 = x88*x97
      x99 = x26*(x14*x67 + x21*x84 + x21*x98)
      x100 = 0.00111534778586481d0*x31
      x101 = x100*x99
      x102 = 4.0d0*x27*x52
      x103 = 0.90856029641607d0*x8
      x104 = 1d0/(x103 + 6.72988144596143d0*x48 + 18.0578d0)
      x105 = 0.01554535d0*log(x103*x104)
      x106 = 0.953184292996937d0*x48
      x107 = x106 + 0.325d0
      x108 = 0.00224786709554261d0*log(x104*x107**2)
      x109 = 1.90636858599387d0*x48
      x110 = x109 + 7.06042d0
      x111 = 0.0524913931697809d0*atan(4.7309269095601136d0/x110)
      x112 = 1d0/(x103 + 3.5529372610885d0*x48 + 12.9352d0)
      x113 = 0.0310907d0*log(x103*x112)
      x114 = x109 + 3.72744d0
      x115 = 0.038783294878113d0*atan(6.1519908197590798d0/x114)
      x116 = x106 + 0.10498d0
      x117 = 0.000969022771154437d0*log(x112*x116**2)
      x118 = x113 + x115 + x117
      x119 = -x105 - x108 - x111 + x118
      x120 = 1.92366105093154d0*x38**1.33333333333333d0 + 1.923661050931
     &54d0*x39** 1.33333333333333d0 - 3.84732210186307d0
      x121 = rho_c**(-4)
      x122 = rho_s**4
      x123 = x121*x122
      x124 = x120*x123
      x125 = x119*x124
      x126 = x109 + 1.13107d0
      x127 = 1d0/(x103 + 1.07811815828004d0*x48 + 13.0045d0)
      x128 = x106 + 0.0047584d0
      x129 = log(x103*x127) + 0.000414033794282063d0*log(x127*x128**2) +
     & 0.317708004743941d0*atan(7.1231089178181177d0/x126)
      x130 = -x123 + 1.0d0
      x131 = x120*x130*x33
      x132 = 0.0974703937105774d0*x131
      x133 = x129*x132
      x134 = x118 - x133
      x135 = -x125 + x134
      x136 = rho_c**(-2.66666666666667d0)
      x137 = 0.0472353356922751d0*mu**8
      x138 = 0.267683468607554d0*x137
      x139 = x136*x138
      x140 = 0.148394183600699d0
      x141 = x25 - 1.0d0
      x142 = x141*x3
      x143 = x140*x142
      x144 = -0.0167302167879722d0*x143*x31
      x145 = rho_c**(-2.0d0)
      x146 = 0.101321183642338d0
      x147 = 0.526297323276602d0*x146
      x148 = x145*x147
      x149 = 4.0d0*mu**6
      x150 = 0.13157433081915d0*x142
      x151 = x15*x17
      x152 = -x150*x151
      x153 = x69*(x38**2.66666666666667d0 + x39**2.66666666666667d0)
      x154 = 0.817704266774463d0*x153
      x155 = 0.558025705063192d0*x23 - 0.352521395009435d0*x8
      x156 = exp(-0.49698248213959023d0*x8)
      x157 = x155*x156*x19
      x158 = x157*x3
      x159 = 10.9027235569928d0*x158
      x160 = 0.00139418473233101d0*x140
      x161 = rho_c**(-1.0d0)*x160
      x162 = 1.55214407110551d0*x13
      x163 = 4.0d0*mu**4
      x164 = 0.508616435555896d0*x23*x42 + 1.0d0
      x165 = x164**(-4)
      x166 = 0.25d0*x165
      x167 = x166*(-x102*(x101 + x3*x62) + x135*x139 + x149*(x135*x148 +
     & x144) + x163 *(x135*x162 + x152 - x161*(-x154 + x159 + x84 + x98)
     &) - x3*x30* x32 - x35*x58)
      x168 = -x113 - x115 - x117
      x169 = x114**(-2)
      x170 = rho_c**(-1.16666666666667d0)*x47
      x171 = x169*x170/(37.8469910464d0*x169 + 1.0d0)
      x172 = 0.30285343213869d0*x10
      x173 = -x172
      x174 = x172*x5
      x175 = x112*(0.592156210181417d0*x170 + x174)
      x176 = 0.90856029641607d0*x7
      x177 = x6*(x173 + x175*x176)
      x178 = (0.000969022771154437d0*x116*x175 - 0.000307885761673592d0*
     &x170)/x116
      x179 = x122/rho_c**5
      x180 = x129*x33
      x181 = x120*x179*x180
      x182 = -x38**0.333333333333333d0 + x39**0.333333333333333d0
      x183 = x130*x180*x182
      x184 = rho_s*x0
      x185 = x105 + x108 + x111 + x168
      x186 = rho_s**5*x185/rho_c**6
      x187 = x179*x185
      x188 = x126**(-2)
      x189 = x127*(0.179686359713341d0*x170 + x174)
      x190 = 0.719040519881222d0*x170*x188/(50.7386806551d0*x188 + 1.0d0
     &) + 1.10064241629821d0*x6*(x173 + x176*x189) + ( 0.000414033794282
     &063d0*x128*x189 - 0.000131550169826529d0*x170)/ x128
      x191 = 0.0758081683534927d0*x171
      x192 = x110**(-2)
      x193 = 0.0342197431724027d0*x177
      x194 = x104*(1.12164690766024d0*x170 + x174)
      x195 = x124*(0.0789023540332771d0*x170*x192/(22.3816694236d0*x192 
     &+ 1.0d0) - x178 - x191 - x193 + 0.0171098715862014d0*x6*(x173 + x1
     &76*x194) + (0.00224786709554261d0*x107*x194 - 0.000714210536071954
     &d0*x170)/ x107)
      x196 = x2 - 1.0d0
      x197 = x196*x25
      x198 = x29*x32
      x199 = x33*(x34 - 1.0d0)
      x200 = x124*x185 + x134
      x201 = x154 - x159 - x84 - x98
      x202 = x20*x31
      x203 = rho_c**(-4.66666666666667d0)
      x204 = 0.17845564573837d0*x30
      x205 = rho_c**(-2.33333333333333d0)
      x206 = x12*x205
      x207 = x10*x5
      x208 = x196*x9*(x145*x18 - 0.045088146776381d0*x202 - 0.0016889954
     &1983155d0* x206 - 0.00627239743302441d0*x207)
      x209 = rho_c**(-3.0d0)
      x210 = x39**(-0.333333333333333d0)
      x211 = x38**(-0.333333333333333d0)
      x212 = x210 - x211
      x213 = x199*x43
      x214 = 2.0d0*x57
      x215 = rho_c**(-3.66666666666667d0)
      x216 = 1.85745434535825d0*x170
      x217 = x207*x50
      x218 = rho_c**(-1.5d0)*x45*x53
      x219 = -x210 + x211
      x220 = rho_s*x219
      x221 = 7.42981738143299d0*rho_c**(-2.16666666666667d0)*x220*x47*x4
     &9
      x222 = rho_c**(-2.33333333333333d0)
      x223 = rho_s*x222
      x224 = mu*x219*x223*x44*x5
      x225 = mu/x55
      x226 = 0.25d0*x183
      x227 = 2.56488140124205d0*x182
      x228 = 4.0d0*x120
      x229 = -x132*x190 + x178 - 0.38988157484231d0*x181 - x184*x226 + x
     &186*x227 - x187*x228 + x191 + x193 + x195
      x230 = rho_c**(-4.0d0)
      x231 = 0.350864882184401d0*x60
      x232 = 24.0d0*x65
      x233 = 1d0/x75
      x234 = x223*x72
      x235 = 0.00706864485168613d0*x16*x68
      x236 = x21*x235
      x237 = 1d0/x89
      x238 = x223*x85
      x239 = rho_c**(-0.333333333333333d0)
      x240 = x239*x69
      x241 = 0.545136177849642d0*x240
      x242 = x38*x78
      x243 = rho_c**(-1.33333333333333d0)*rho_s
      x244 = x243*x69
      x245 = 2.18054471139857d0*x244
      x246 = x39*x95
      x247 = x246*x87
      x248 = 0.034943218589452d0*rho_s*x136*x20
      x249 = 0.164800693786261d0*x5
      x250 = 0.034943218589452d0*x202
      x251 = 0.164800693786261d0*x207
      x252 = x250*x73 + x251*x75
      x253 = x21*x70/x77**2
      x254 = x250*x86 + x251*x89
      x255 = x21*x88/x94**2
      x256 = -x233*x236*x79*(x10 - x234) - x236*x237*x96*(x10 + x238) - 
     &x241*x83 - x241*x87*x97 + x242*x245*x82 - x245*x247*x92 - x253*x71
     &*x82*( -x234*x249*x75 - x248*x72*x73 + x252) - x255*x92*x93*(x238*
     &x249* x89 + x248*x85*x86 + x254)
      x257 = x100*x26
      x258 = x1*x141
      x259 = 0.0334604335759443d0*x140
      x260 = x145*x17
      x261 = x3*x9*(-0.0225440733881905d0*x202 - 0.000844497709915773d0*
     &x206 + 0.113935126864452d0*x207*x24 - 0.0031361987165122d0*x207 + 
     &0.0047892375d0*x260)
      x262 = 0.263148661638301d0*x17
      x263 = x156*x3
      x264 = 21.8054471139857d0*x157
      x265 = -2.18054471139857d0*x38**1.66666666666667d0 + 2.18054471139
     &857d0*x39** 1.66666666666667d0
      x266 = x160*x31
      x267 = rho_s**3
      x268 = x121*x267
      x269 = 0.38988157484231d0*x120*x180
      x270 = x119*x227
      x271 = x119*x228
      x272 = 7.42981738143299d0*x170
      x273 = x267/rho_c**3
      x274 = x123*x270 + x226 + x269*x273 - x271*x273
      x275 = rho_s*x141
      x276 = x22*x235
      x277 = 2.18054471139857d0*x240
      x278 = x80*(x81 - 1.0d0)
      x279 = x91 - 1.0d0
      x280 = -x233*x242*x276 + x237*x246*x276 + x242*x277*x278 - x247*x2
     &77*x279 + x252*x253*x278*x38 - x254*x255*x279*x39
      E = -rho_c*(x125 + x133 + x167 + x168)
      d1E(1) = -0.25d0*rho_c*(16.0d0*x120*x187 + 0.38988157484231d0*x131
     &*x190 + x165*( -x1*x203*x204 + x102*(-x1*x230*x231 + 0.00185891297
     &644135d0*x136* x99 - 0.175432441092201d0*x145*x196*x60 + 0.0049018
     &0588786583d0* x197*x205 + x208*x59*x61 + x257*(-2.41662179562687d0
     &*rho_c**( -0.333333333333333d0)*x3*x63*x64 - x1*x145*x232 - x14*x1
     &6*x64*x66 *(-0.550321208149104d0*x202 + 1.49912448908652d0*x207) +
     & x256 - x67)) - 0.148713038115308d0*x136*x197*x29 - 0.713822582953
     &479d0* x137*x200*x215 + x139*x229 + x149*(0.0278836946466203d0*x13
     &6*x143 - 1.0525946465532d0*x146*x200*x209 + x148*x229 - x203*x258*
     &x259 - x259*x261*x31) + x163*(-0.00232364122055169d0*x140*x145*x20
     &1 + x150*x260 - 0.263148661638301d0*x151*x261 + x162*x229 - 2.0695
     &2542814068d0*x200*x206 - x230*x258*x262 + x266*(-x1*x222* x264 + 0
     &.545136177849642d0*x153*x239 - 1.80615420514536d0*x155* x22*x263*x
     &4 - 7.26848237132856d0*x158*x239 - 10.9027235569928d0* x19*x21*x26
     &3*(-0.372017136708795d0*x202 + 0.117507131669812d0* x207) + x244*x
     &265 + x256)) + x184*x212*x213*x214 + 0.00116228665296198d0*x197*x2
     &09*x28 + x198*x208 + x213*x225*( 54.281640924821d0*rho_c**(-2.5d0)
     &*x220*x42*x54 - x216 - 9.02445731974005d0*x217 - 13.5704102312052d
     &0*x218 + x221 + 36.0978292789602d0*x224 + x56*(x216 + 4.1775723570
     &5837d0*x217 - x221 - 16.7102894282335d0*x224))) - 0.30323267341397
     &1d0*x171 - 0.136878972689611d0*x177 - 4.0d0*x178 + 1.5595262993692
     &4d0*x181 - 10.2595256049682d0*x182*x186 + x183*x184 - 4.0d0*x195 +
     & 1.35631049481572d0*x202*x42*(x102*(-x101 + x196*x62) + x139*x200 
     &+ x149*(x144 + x148*x200) + x163*(x152 + x161*x201 + x162*x200) + 
     &x197*x198 + x199*x58)/x164**5) + x135 - x167
      d1E(2) = rho_c*(-x166*(rho_s*x204*x215 + x102*(rho_s*x209*x231 + x
     &257*(rho_s*x15* x232 + x280)) + x138*x215*x274 + x149*(x147*x209*x
     &274 + x215*x259 *x275) + x163*(1.55214407110551d0*x206*x274 + x209
     &*x262*x275 + x266*(-x240*x265 + x243*x264 + x280)) - x214*x219*x35
     &*x36*x43 + x225*x35*x40*(x212*x56*(16.7102894282335d0*x217 + x272)
     & + x219*( 36.0978292789602d0*x217 + 54.281640924821d0*x218 + x272)
     &)) + x179 *x270 + x226*x36 + x268*x269 - x268*x271)
      end subroutine


C*****************************************************************************
      pure subroutine D2ESRC_SPIN_VWN5_ERF(rho_c, rho_s, mu, E, d1E, d2E
     &)
C*****************************************************************************
C   Implemented by E.R. Kjellgren.
C
C   Subroutine generated using Sympy 1.3
C   Generated: March 25, 2019
C*****************************************************************************
      implicit none
      real*8, intent(in) :: rho_c, rho_s, mu
      real*8, intent(out) :: E, d1E(9), d2E(45)
      real*8 :: x0, x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12, x
     &13, x14, x15, x16, x17, x18, x19, x20, x21, x22, x23, x24, x25, x2
     &6, x27, x28, x29, x30, x31, x32, x33, x34, x35, x36, x37, x38, x39
     &, x40, x41, x42, x43, x44, x45, x46, x47, x48, x49, x50, x51, x52,
     & x53, x54, x55, x56, x57, x58, x59, x60, x61, x62, x63, x64, x65, 
     &x66, x67, x68, x69, x70, x71, x72, x73, x74, x75, x76, x77, x78, x
     &79, x80, x81, x82, x83, x84, x85, x86, x87, x88, x89, x90, x91, x9
     &2, x93, x94, x95, x96, x97, x98, x99, x100, x101, x102, x103, x104
     &, x105, x106, x107, x108, x109, x110, x111, x112, x113, x114, x115
     &, x116, x117, x118, x119, x120, x121, x122, x123, x124, x125, x126
     &, x127, x128, x129, x130, x131, x132, x133, x134, x135, x136, x137
     &, x138, x139, x140, x141, x142, x143, x144, x145, x146, x147, x148
     &, x149, x150, x151, x152, x153, x154, x155, x156, x157, x158, x159
     &, x160, x161, x162, x163, x164, x165, x166, x167, x168, x169, x170
     &, x171, x172, x173, x174, x175, x176, x177, x178, x179, x180, x181
     &, x182, x183, x184, x185, x186, x187, x188, x189, x190, x191, x192
     &, x193, x194, x195, x196, x197, x198, x199, x200, x201, x202, x203
     &, x204, x205, x206, x207, x208, x209, x210, x211, x212, x213, x214
     &, x215, x216, x217, x218, x219, x220, x221, x222, x223, x224, x225
     &, x226, x227, x228, x229, x230, x231, x232, x233, x234, x235, x236
     &, x237, x238, x239, x240, x241, x242, x243, x244, x245, x246, x247
     &, x248, x249, x250, x251, x252, x253, x254, x255, x256, x257, x258
     &, x259, x260, x261, x262, x263, x264, x265, x266, x267, x268, x269
     &, x270, x271, x272, x273, x274, x275, x276, x277, x278, x279, x280
     &, x281, x282, x283, x284, x285, x286, x287, x288, x289, x290, x291
     &, x292, x293, x294, x295, x296, x297, x298, x299, x300, x301, x302
     &, x303, x304, x305, x306, x307, x308, x309, x310, x311, x312, x313
     &, x314, x315, x316, x317, x318, x319, x320, x321, x322, x323, x324
     &, x325, x326, x327, x328, x329, x330, x331, x332, x333, x334, x335
     &, x336, x337, x338, x339, x340, x341, x342, x343, x344, x345, x346
     &, x347, x348, x349, x350, x351, x352, x353, x354, x355, x356, x357
     &, x358, x359, x360, x361, x362, x363, x364, x365, x366, x367, x368
     &, x369, x370, x371, x372, x373, x374, x375, x376, x377, x378, x379
     &, x380, x381, x382, x383, x384, x385, x386, x387, x388, x389, x390
     &, x391, x392, x393, x394, x395, x396, x397, x398, x399, x400, x401
     &, x402, x403, x404, x405, x406, x407, x408, x409, x410, x411, x412
     &, x413, x414, x415, x416, x417, x418, x419, x420, x421, x422, x423
     &, x424, x425, x426, x427, x428, x429, x430, x431, x432, x433, x434
     &, x435, x436, x437, x438, x439, x440, x441, x442, x443, x444, x445
     &, x446, x447, x448, x449, x450, x451, x452, x453, x454, x455, x456
     &, x457, x458, x459, x460, x461, x462, x463, x464, x465, x466, x467
     &, x468, x469, x470, x471, x472, x473, x474, x475, x476, x477, x478
     &, x479, x480, x481, x482, x483, x484, x485, x486, x487, x488, x489
     &, x490, x491, x492, x493, x494, x495, x496, x497, x498, x499, x500
     &, x501, x502, x503, x504, x505, x506, x507, x508, x509, x510, x511
     &, x512, x513, x514, x515, x516, x517, x518, x519, x520, x521, x522
     &, x523, x524, x525, x526, x527, x528, x529, x530, x531, x532, x533
     &, x534, x535, x536, x537, x538, x539, x540, x541, x542, x543, x544
     &, x545, x546, x547, x548, x549, x550, x551, x552, x553, x554, x555
     &, x556, x557, x558, x559, x560, x561, x562, x563, x564, x565, x566
     &, x567, x568, x569, x570, x571, x572, x573, x574, x575, x576, x577
     &, x578, x579, x580, x581, x582, x583, x584, x585, x586, x587, x588
     &, x589, x590, x591, x592, x593, x594, x595, x596, x597, x598, x599
     &, x600, x601, x602, x603, x604, x605, x606, x607, x608, x609, x610
     &, x611, x612, x613, x614, x615, x616, x617, x618, x619, x620, x621
     &, x622, x623, x624, x625, x626, x627, x628, x629, x630, x631, x632
     &, x633, x634, x635, x636, x637, x638, x639, x640, x641, x642, x643
     &, x644, x645, x646, x647, x648, x649, x650, x651, x652, x653, x654
     &, x655, x656, x657, x658, x659, x660, x661, x662, x663, x664, x665
     &, x666, x667, x668, x669, x670, x671, x672, x673, x674, x675, x676
     &, x677, x678, x679, x680, x681, x682, x683, x684, x685, x686, x687
     &, x688, x689, x690, x691, x692, x693, x694, x695, x696, x697, x698
     &, x699, x700, x701, x702, x703, x704, x705, x706, x707, x708, x709
     &, x710, x711, x712, x713, x714, x715, x716, x717, x718, x719, x720
     &, x721, x722, x723, x724, x725, x726, x727, x728, x729, x730, x731
     &, x732, x733, x734, x735, x736, x737, x738, x739, x740, x741, x742
     &, x743, x744, x745, x746, x747, x748, x749, x750, x751, x752, x753
     &, x754, x755, x756, x757, x758, x759, x760, x761, x762, x763, x764
     &, x765, x766, x767, x768, x769, x770, x771, x772, x773, x774, x775
     &, x776, x777, x778, x779, x780, x781, x782, x783, x784, x785, x786
     &, x787, x788, x789, x790, x791, x792, x793, x794, x795, x796, x797
     &, x798, x799, x800, x801, x802, x803, x804, x805, x806, x807, x808
     &
      E = 0.0d0
      d1E(:) = 0.0d0
      d2E(:) = 0.0d0
      x0 = rho_c**(-2)
      x1 = rho_s**2
      x2 = x0*x1
      x3 = -x2 + 1.0d0
      x4 = 1.46459188756152d0
      x5 = 1d0/x4
      x6 = rho_c**0.333333333333333d0
      x7 = 1d0/x6
      x8 = x5*x7
      x9 = exp(-0.6836107611867116d0*x8)
      x10 = rho_c**(-1.33333333333333d0)
      x11 = 4.60115111447049d0
      x12 = 1d0/x11
      x13 = x10*x12
      x14 = rho_c**1.0d0
      x15 = 1d0/x14
      x16 = 3.14159265358979d0
      x17 = 1d0/x16
      x18 = 0.009578475d0*x17
      x19 = 2.14502939711103d0
      x20 = 1d0/x19
      x21 = rho_c**0.666666666666667d0
      x22 = 1d0/x21
      x23 = x20*x22
      x24 = 0.00126674656487366d0*x13 - x15*x18 + 0.0676322201645715d0*x
     &23 + 0.0188171922990732d0*x8 + 1.0d0
      x25 = x24*x9
      x26 = x25*x3
      x27 = 0.0837224526465878d0
      x28 = 1.41421356237310d0
      x29 = mu**5*x28
      x30 = x27*x29
      x31 = rho_c**(-1.66666666666667d0)
      x32 = 0.0892278228691848d0*x31
      x33 = x30*x32
      x34 = 0.101321183642338d0
      x35 = 0.693147180559945d0
      x36 = x34*(-x35 + 1.0d0)
      x37 = 1d0/rho_c
      x38 = rho_s*x37
      x39 = x38 + 1.0d0
      x40 = x39**0.666666666666667d0
      x41 = -x38 + 1.0d0
      x42 = x41**0.666666666666667d0
      x43 = x40 + x42
      x44 = x43**3
      x45 = mu**2
      x46 = x43**2
      x47 = 1d0/x46
      x48 = x45*x47
      x49 = x48*x8
      x50 = 0.826307487110758d0
      x51 = rho_c**(-0.166666666666667d0)*x50
      x52 = 1d0/x43
      x53 = mu*x52
      x54 = 11.1447260721495d0*x51*x53 + 1.0d0
      x55 = 12.5327170711751d0*x49 + x54
      x56 = 1d0/x55
      x57 = mu**3
      x58 = 0.564189583547756d0
      x59 = 1d0/x44
      x60 = x58*x59
      x61 = 27.1408204624105d0*rho_c**(-0.5d0)*x57*x60 + 27.073371959220
     &1d0*x49 + x54
      x62 = x56*x61
      x63 = log(x62)
      x64 = x44*x63
      x65 = 0.179587122125167d0
      x66 = x25*x65
      x67 = 0.175432441092201d0*x15
      x68 = x66*x67
      x69 = 0.825481812223657d0*x23 - 4.49737346725955d0*x8
      x70 = exp(-0.28165369188898165d0*x8)
      x71 = x16*x70
      x72 = x3*x71
      x73 = 12.0d0*x72
      x74 = x69*x73
      x75 = 2.28942848510666d0
      x76 = x11*x75
      x77 = 0.817704266774463d0*x76
      x78 = x39**2
      x79 = 1d0/x39
      x80 = x79**0.666666666666667d0
      x81 = 0.0524148278841779d0*x23
      x82 = x79**0.333333333333333d0
      x83 = x8*x82
      x84 = x80*x81 + 0.494402081358784d0*x83 + 1.0d0
      x85 = 1d0/x84
      x86 = 1d0/x80
      x87 = 0.0259335011650457d0*x83
      x88 = x86*(-x87 + 1.0d0)
      x89 = x85*x88
      x90 = x78*x89
      x91 = x77*x90
      x92 = x41**2
      x93 = 1d0/x41
      x94 = x93**0.333333333333333d0
      x95 = x8*x94
      x96 = 0.0259335011650457d0*x95
      x97 = -x96 + 1.0d0
      x98 = x93**0.666666666666667d0
      x99 = x81*x98 + 0.494402081358784d0*x95 + 1.0d0
      x100 = 1d0/x99
      x101 = 1d0/x98
      x102 = x100*x101
      x103 = x102*x97
      x104 = x103*x92
      x105 = x104*x77
      x106 = x27*(x105*x21 + x14*x74 + x21*x91)
      x107 = 0.00111534778586481d0*x31
      x108 = x106*x107
      x109 = x28*x57
      x110 = 4.0d0*x109
      x111 = 0.90856029641607d0*x8
      x112 = x111 + 6.72988144596143d0*x51 + 18.0578d0
      x113 = 1d0/x112
      x114 = 0.01554535d0*log(x111*x113)
      x115 = 0.953184292996937d0*x51
      x116 = x115 + 0.325d0
      x117 = x116**2
      x118 = x113*x117
      x119 = 0.00224786709554261d0*log(x118)
      x120 = 1.90636858599387d0*x51
      x121 = x120 + 7.06042d0
      x122 = 0.0524913931697809d0*atan(4.7309269095601136d0/x121)
      x123 = x111 + 3.5529372610885d0*x51 + 12.9352d0
      x124 = 1d0/x123
      x125 = 0.0310907d0*log(x111*x124)
      x126 = x120 + 3.72744d0
      x127 = 0.038783294878113d0*atan(6.1519908197590798d0/x126)
      x128 = x115 + 0.10498d0
      x129 = x128**2
      x130 = x124*x129
      x131 = 0.000969022771154437d0*log(x130)
      x132 = x125 + x127 + x131
      x133 = -x114 - x119 - x122 + x132
      x134 = 1.92366105093154d0*x39**1.33333333333333d0 + 1.923661050931
     &54d0*x41** 1.33333333333333d0 - 3.84732210186307d0
      x135 = rho_c**(-4)
      x136 = rho_s**4
      x137 = x135*x136
      x138 = x134*x137
      x139 = x133*x138
      x140 = x120 + 1.13107d0
      x141 = x111 + 1.07811815828004d0*x51 + 13.0045d0
      x142 = 1d0/x141
      x143 = x115 + 0.0047584d0
      x144 = x143**2
      x145 = x142*x144
      x146 = 0.000414033794282063d0*log(x145) + log(x111*x142) + 0.31770
     &8004743941d0* atan(7.1231089178181177d0/x140)
      x147 = x34*(-x137 + 1.0d0)
      x148 = x134*x147
      x149 = 0.0974703937105774d0*x148
      x150 = x146*x149
      x151 = x132 - x150
      x152 = -x139 + x151
      x153 = rho_c**(-2.66666666666667d0)
      x154 = 0.0472353356922751d0*mu**8
      x155 = 0.267683468607554d0*x154
      x156 = x153*x155
      x157 = 0.148394183600699d0
      x158 = x25 - 1.0d0
      x159 = x157*x158
      x160 = x159*x3
      x161 = -0.0167302167879722d0*x160*x31
      x162 = rho_c**(-2.0d0)
      x163 = 0.101321183642338d0
      x164 = 0.526297323276602d0*x163
      x165 = x162*x164
      x166 = mu**6
      x167 = 4.0d0*x166
      x168 = x158*x3
      x169 = 0.13157433081915d0*x168
      x170 = x15*x17
      x171 = -x169*x170
      x172 = x76*(x39**2.66666666666667d0 + x41**2.66666666666667d0)
      x173 = 0.817704266774463d0*x172
      x174 = exp(-0.49698248213959023d0*x8)
      x175 = x174*x19
      x176 = 0.558025705063192d0*x23 - 0.352521395009435d0*x8
      x177 = x176*x3
      x178 = x175*x177
      x179 = 10.9027235569928d0*x178
      x180 = 0.00139418473233101d0*x157
      x181 = rho_c**(-1.0d0)*x180
      x182 = 1.55214407110551d0*x13
      x183 = mu**4
      x184 = 4.0d0*x183
      x185 = -x110*(x108 + x3*x68) + x152*x156 + x167*(x152*x165 + x161)
     & + x184*(x152 *x182 + x171 - x181*(x105 - x173 + x179 + x91)) - x2
     &6*x33 - x36* x64
      x186 = 0.508616435555896d0*x23*x45 + 1.0d0
      x187 = x186**(-4)
      x188 = 0.25d0*x187
      x189 = x185*x188
      x190 = -x125 - x127 - x131
      x191 = rho_c**(-1.16666666666667d0)*x50
      x192 = x126**(-2)
      x193 = 37.8469910464d0*x192 + 1.0d0
      x194 = 1d0/x193
      x195 = x192*x194
      x196 = x191*x195
      x197 = 0.30285343213869d0*x10
      x198 = -x197
      x199 = x197*x5
      x200 = 0.592156210181417d0*x191 + x199
      x201 = x124*x200
      x202 = 0.90856029641607d0*x7
      x203 = x201*x202
      x204 = x198 + x203
      x205 = x204*x6
      x206 = 1d0/x128
      x207 = 0.000307885761673592d0*x191
      x208 = x128*x201
      x209 = 0.000969022771154437d0*x208
      x210 = x206*(-x207 + x209)
      x211 = rho_c**(-5)
      x212 = x136*x211
      x213 = x212*x34
      x214 = x134*x146
      x215 = x213*x214
      x216 = x41**0.333333333333333d0
      x217 = x39**0.333333333333333d0
      x218 = x216 - x217
      x219 = x146*x218
      x220 = rho_s*x0
      x221 = x147*x220
      x222 = x219*x221
      x223 = rho_c**(-6)
      x224 = x114 + x119 + x122 + x190
      x225 = rho_s**5
      x226 = x224*x225
      x227 = x223*x226
      x228 = 10.2595256049682d0*x218
      x229 = x212*x224
      x230 = x134*x229
      x231 = x140**(-2)
      x232 = 50.7386806551d0*x231 + 1.0d0
      x233 = 1d0/x232
      x234 = x231*x233
      x235 = 0.719040519881222d0*x191*x234
      x236 = 0.179686359713341d0*x191 + x199
      x237 = x142*x236
      x238 = x202*x237
      x239 = x198 + x238
      x240 = 1.10064241629821d0*x6
      x241 = x239*x240
      x242 = 1d0/x143
      x243 = 0.000131550169826529d0*x191
      x244 = x143*x237
      x245 = 0.000414033794282063d0*x244
      x246 = x235 + x241 + x242*(-x243 + x245)
      x247 = 0.38988157484231d0*x246
      x248 = 1d0/x116
      x249 = 0.000714210536071954d0*x191
      x250 = 1.12164690766024d0*x191 + x199
      x251 = x113*x250
      x252 = x116*x251
      x253 = 0.00224786709554261d0*x252
      x254 = x202*x251
      x255 = x198 + x254
      x256 = 0.0171098715862014d0*x6
      x257 = x255*x256
      x258 = 0.0342197431724027d0*x205
      x259 = x121**(-2)
      x260 = 22.3816694236d0*x259 + 1.0d0
      x261 = 1d0/x260
      x262 = x259*x261
      x263 = 0.0758081683534927d0*x196
      x264 = 0.0789023540332771d0*x191*x262 - x263
      x265 = -x210 + x248*(-x249 + x253) + x257 - x258 + x264
      x266 = x138*x265
      x267 = x2 - 1.0d0
      x268 = x25*x267
      x269 = x268*x27
      x270 = x269*x29
      x271 = x34*(x35 - 1.0d0)
      x272 = x138*x224 + x151
      x273 = -x105 + x173 - x179 - x91
      x274 = x110*(-x108 + x267*x68) + x156*x272 + x167*(x161 + x165*x27
     &2) + x184*( x171 + x181*x273 + x182*x272) + x270*x32 + x271*x64
      x275 = x186**(-5)
      x276 = x20*x31
      x277 = x275*x276*x45
      x278 = rho_c**(-4.66666666666667d0)
      x279 = x1*x278
      x280 = x25*x30
      x281 = 0.17845564573837d0*x280
      x282 = -x279*x281
      x283 = rho_c**(-2.33333333333333d0)
      x284 = x12*x283
      x285 = 0.00168899541983155d0*x284
      x286 = x162*x18
      x287 = 0.045088146776381d0*x276
      x288 = x10*x5
      x289 = 0.00627239743302441d0*x288
      x290 = -x285 + x286 - x287 - x289
      x291 = x267*x9
      x292 = x290*x291
      x293 = x268*x29
      x294 = rho_c**(-3.0d0)
      x295 = 0.0571643564037363d0
      x296 = x294*x295
      x297 = 0.0203323666368788d0*x296
      x298 = x153*x27
      x299 = 0.148713038115308d0*x298
      x300 = x271*x46
      x301 = 2.0d0*x63
      x302 = x300*x301
      x303 = x41**(-0.333333333333333d0)
      x304 = x39**(-0.333333333333333d0)
      x305 = x303 - x304
      x306 = x220*x305
      x307 = x154*x272
      x308 = rho_c**(-3.66666666666667d0)
      x309 = 0.713822582953479d0*x308
      x310 = x288*x53
      x311 = 4.17757235705837d0*x310
      x312 = mu*x5
      x313 = x312*x47
      x314 = -x303 + x304
      x315 = rho_c**(-2.33333333333333d0)
      x316 = rho_s*x315
      x317 = x314*x316
      x318 = x313*x317
      x319 = 16.7102894282335d0*x318
      x320 = 1.85745434535825d0*x191
      x321 = rho_c**(-2.16666666666667d0)*x50
      x322 = x321*x52
      x323 = rho_s*x314
      x324 = x322*x323
      x325 = 7.42981738143299d0*x324
      x326 = x320 - x325
      x327 = x311 - x319 + x326
      x328 = 9.02445731974005d0*x310
      x329 = x48*x58
      x330 = rho_c**(-1.5d0)*x329
      x331 = 13.5704102312052d0*x330
      x332 = 36.0978292789602d0*x318
      x333 = rho_c**(-2.5d0)
      x334 = x45*x60
      x335 = x333*x334
      x336 = 54.281640924821d0*x323*x335
      x337 = -x320 + x325
      x338 = -x328 - x331 + x332 + x336 + x337
      x339 = x327*x62 + x338
      x340 = 1d0/x61
      x341 = mu*x340
      x342 = x339*x341
      x343 = 2.56488140124205d0*x218
      x344 = 0.25d0*x147
      x345 = x219*x344
      x346 = -0.38988157484231d0*x215 - x220*x345 + x263
      x347 = -x149*x246 + x210 + x227*x343 - 4.0d0*x230 + x258 + x266 + 
     &x346
      x348 = -0.550321208149104d0*x276 + 1.49912448908652d0*x288
      x349 = x14*x73
      x350 = x348*x349
      x351 = x69*x71
      x352 = 24.0d0*x351
      x353 = x1*x352
      x354 = x162*x353
      x355 = x3*x69*x70
      x356 = 2.14502939711103d0
      x357 = rho_c**(-0.333333333333333d0)*x356
      x358 = 1.12661476755593d0*x355*x357
      x359 = 0.00706864485168613d0*x21
      x360 = x316*x79
      x361 = x10 - x360
      x362 = x16*x75
      x363 = 1d0/x82
      x364 = x363*x85
      x365 = x362*x364
      x366 = x361*x365*x78
      x367 = x359*x366
      x368 = 1d0/x94
      x369 = x100*x368
      x370 = x362*x369
      x371 = x370*x92
      x372 = x316*x93
      x373 = x10 + x372
      x374 = x359*x373
      x375 = x371*x374
      x376 = rho_c**(-0.333333333333333d0)
      x377 = x376*x76
      x378 = 0.545136177849642d0*x377
      x379 = x378*x90
      x380 = x104*x378
      x381 = x39*x89
      x382 = rho_c**(-1.33333333333333d0)
      x383 = x382*x76
      x384 = rho_s*x383
      x385 = 2.18054471139857d0*x384
      x386 = x381*x385
      x387 = x103*x41
      x388 = x385*x387
      x389 = x84**(-2)
      x390 = x389*x88
      x391 = x21*x77
      x392 = x390*x391
      x393 = x153*x20
      x394 = 0.034943218589452d0*x393
      x395 = rho_s*x79
      x396 = x395*x80
      x397 = 0.164800693786261d0*x5
      x398 = x360*x82
      x399 = 0.034943218589452d0*x276
      x400 = 0.164800693786261d0*x288
      x401 = x399*x80 + x400*x82
      x402 = -x394*x396 - x397*x398 + x401
      x403 = x402*x78
      x404 = x392*x403
      x405 = rho_s*x93
      x406 = x405*x98
      x407 = x372*x94
      x408 = x399*x98 + x400*x94
      x409 = x394*x406 + x397*x407 + x408
      x410 = x99**(-2)
      x411 = x101*x97
      x412 = x410*x411
      x413 = x391*x412
      x414 = x413*x92
      x415 = x409*x414
      x416 = -x367 - x375 - x379 - x380 + x386 - x388 - x404 - x415
      x417 = -x350 - x354 - x358 + x416 - x74
      x418 = x107*x27
      x419 = 0.122619224952946d0
      x420 = x283*x419
      x421 = 0.0399758348639607d0*x420
      x422 = x65*x67
      x423 = 0.175432441092201d0*x162*x66
      x424 = 0.00185891297644135d0*x153
      x425 = 0.350864882184401d0*x66
      x426 = rho_c**(-4.0d0)
      x427 = x1*x426
      x428 = x106*x424 - x425*x427
      x429 = 0.0334604335759443d0*x159
      x430 = x279*x429
      x431 = x153*x157
      x432 = 0.0278836946466203d0*x168*x431
      x433 = x162*x17
      x434 = -0.0225440733881905d0*x276 - 0.000844497709915773d0*x284 - 
     &0.0031361987165122d0*x288 + 0.0047892375d0*x433
      x435 = 0.113935126864452d0*x24*x288 + x434
      x436 = x3*x9
      x437 = x435*x436
      x438 = 0.0334604335759443d0*x157*x31
      x439 = x437*x438
      x440 = x163*x272
      x441 = 1.0525946465532d0*x294*x440
      x442 = x165*x347
      x443 = 0.263148661638301d0*x158
      x444 = x1*x17
      x445 = x426*x443*x444
      x446 = x169*x433
      x447 = 0.263148661638301d0*x170
      x448 = x437*x447
      x449 = x157*x273
      x450 = 0.00232364122055169d0*x162*x449
      x451 = 2.06952542814068d0*x272*x284
      x452 = -0.372017136708795d0*x276 + 0.117507131669812d0*x288
      x453 = x175*x452
      x454 = x3*x453
      x455 = 10.9027235569928d0*x21
      x456 = x175*x176
      x457 = 21.8054471139857d0*x456
      x458 = x1*x315
      x459 = x174*x177
      x460 = x22*x4
      x461 = x41**1.66666666666667d0
      x462 = x39**1.66666666666667d0
      x463 = x76*(x461 - x462)
      x464 = rho_s*x382
      x465 = 0.545136177849642d0*x172*x376 - 7.26848237132856d0*x178*x37
     &6 + x416 - x454*x455 - x457*x458 - 1.80615420514536d0*x459*x460 + 
     &2.18054471139857d0*x463*x464
      x466 = x180*x31
      x467 = x465*x466
      x468 = x182*x347
      x469 = x110*(-x267*x423 + x268*x421 + x292*x422 + x417*x418 + x428
     &) + x156*x347 + x167*(-x430 + x432 - x439 - x441 + x442) + x184*(-
     &x445 + x446 - x448 - x450 - x451 + x467 + x468) + x282 + x292*x33 
     &+ x293*x297 - x293*x299 + x300*x342 + x302*x306 - x307*x309
      x470 = rho_s**3
      x471 = x135*x470
      x472 = x214*x34
      x473 = 0.38988157484231d0*x472
      x474 = x133*x343
      x475 = x133*x134
      x476 = 4.0d0*x475
      x477 = rho_s*x308
      x478 = x281*x477
      x479 = x314*x37
      x480 = x36*x46
      x481 = x301*x480
      x482 = x341*x43
      x483 = 7.42981738143299d0*x191
      x484 = 36.0978292789602d0*x310 + 54.281640924821d0*x330 + x483
      x485 = x314*x484
      x486 = 16.7102894282335d0*x310 + x483
      x487 = x486*x62
      x488 = x36*(x305*x487 + x485)
      x489 = rho_c**(-3)
      x490 = x470*x489
      x491 = x473*x490
      x492 = x137*x474 + x345 - x476*x490 + x491
      x493 = x155*x308
      x494 = x308*x429
      x495 = rho_s*x494
      x496 = x164*x294
      x497 = x294*x425
      x498 = rho_s*x497
      x499 = x15*x352
      x500 = x369*x41
      x501 = x22*x362
      x502 = 0.00706864485168613d0*x501
      x503 = x364*x39
      x504 = x500*x502 - x502*x503
      x505 = rho_s*x499 + x504
      x506 = 2.18054471139857d0*x377
      x507 = x39*x506
      x508 = x86*(x87 - 1.0d0)
      x509 = x96 - 1.0d0
      x510 = x41*x506
      x511 = x39*x391
      x512 = x391*x41
      x513 = -x101*x408*x410*x509*x512 - x102*x509*x510 + x389*x401*x508
     &*x511 + x507* x508*x85
      x514 = x17*x294
      x515 = x443*x514
      x516 = rho_s*x515
      x517 = 1.55214407110551d0*x284
      x518 = 2.18054471139857d0*x376
      x519 = x382*x457
      x520 = rho_s*x519 + x504
      x521 = -x188*(x110*(x418*(x505 + x513) + x498) + x167*(x492*x496 +
     & x495) + x184 *(x466*(-x463*x518 + x513 + x520) + x492*x517 + x516
     &) + x478 - x479*x481 + x482*x488 + x492*x493) + x212*x474 + x345*x
     &37 + x471* x473 - x471*x476
      x522 = x197 - x203
      x523 = 0.000615771523347183d0*x191
      x524 = x223*x225
      x525 = x133*x218
      x526 = 8.0d0*x212
      x527 = x235 - x240*(x197 - x238) - x242*(x243 - x245)
      x528 = x206*(x207 - x209)
      x529 = 0.0342197431724027d0*x6
      x530 = x522*x529
      x531 = x138*(-x248*(x249 - x253) - x256*(x197 - x254) + x264 + x52
     &8 + x530)
      x532 = x436*(x285 - x286 + x287 + x289)
      x533 = x26*x29
      x534 = x0*x481
      x535 = x315*x5
      x536 = x535/(x126**5*x193**2)
      x537 = x194*x535/x126**3
      x538 = rho_c**(-2.16666666666667d0)*x50
      x539 = x195*x538
      x540 = rho_c**(-0.666666666666667d0)
      x541 = x204*x540
      x542 = 1d0/x129
      x543 = -0.317728097665645d0*x191
      x544 = x208 + x543
      x545 = x542*x544
      x546 = x201*x206*x544
      x547 = x219*x34
      x548 = rho_c**(-7)
      x549 = x225*x548
      x550 = x136*x223
      x551 = x472*x550
      x552 = 0.40380457618492d0*x283
      x553 = 0.60570686427738d0*x10
      x554 = x5*x552
      x555 = 0.69084891187832d0*x538 + x554
      x556 = 0.60570686427738d0*x288
      x557 = x200*(1.18431242036283d0*x191 + x556)/x123**2
      x558 = -x124*x202*x555 - x201*x553 + x202*x557 + x552
      x559 = x39**(-0.666666666666667d0)
      x560 = 0.854960467080683d0*x38
      x561 = x559*x560
      x562 = x41**(-0.666666666666667d0)
      x563 = x560*x562
      x564 = -5.1297628024841d0*x216 + 5.1297628024841d0*x217 + x561 + x
     &563
      x565 = x146*x147
      x566 = rho_s*x489
      x567 = x564*x565*x566
      x568 = x226*x548
      x569 = x218*x568
      x570 = x134*x246
      x571 = x213*x570
      x572 = x134*x224
      x573 = x550*x572
      x574 = x218*x221*x246
      x575 = x542*(0.000359200055285857d0*x128*x538 + 0.0009690227711544
     &37d0*x129* x557 - 0.000969022771154437d0*x130*x555 - x208*x523 + 4
     &.89119786774443d-5*x535)
      x576 = x564*x568
      x577 = x218*x265*x524
      x578 = x134*x265
      x579 = 1d0/x144
      x580 = x244 + x543
      x581 = 0.209634086332231d0*x538 + x554
      x582 = x236*(0.359372719426682d0*x191 + x556)/x141**2
      x583 = -0.838880606528093d0*x234*x538 - x237*x241 - 0.000414033794
     &282063d0*x237 *x242*x580 + 0.366880805432736d0*x239*x540 + x240*(-
     &x142*x202* x581 + x202*x582 - x237*x553 + x552) + x243*x579*x580 +
     & x579*( 0.000153475198130951d0*x143*x538 + 0.000414033794282063d0*
     &x144* x582 - 0.000414033794282063d0*x145*x581 - 0.0002631003396530
     &58d0* x191*x244 + 2.08985926032878d-5*x535) + 0.456918753052755d0*
     &x233* x535/x140**3 - 23.1834546964702d0*x535/(x140**5*x232**2)
      x584 = 1.82319440383792d0*x536
      x585 = 0.0481727702369445d0*x537
      x586 = 0.0884428630790749d0*x539
      x587 = 0.0114065810574676d0*x541
      x588 = x207*x545
      x589 = 1d0/x117
      x590 = x252 + x543
      x591 = x201*x258
      x592 = 0.000969022771154437d0*x546
      x593 = x529*x558
      x594 = 1.30858805893694d0*x538 + x554
      x595 = x250*(2.24329381532048d0*x191 + x556)/x112**2
      x596 = x138*(-0.00224786709554261d0*x248*x251*x590 + x249*x589*x59
     &0 - x251*x257 + 0.00570329052873379d0*x255*x540 + x256*(-x113*x202
     &*x594 + x202* x595 - x251*x553 + x552) - 0.0920527463721566d0*x262
     &*x538 - x575 + x584 - x585 + x586 - x587 - x588 + x589*(0.00083324
     &5625417279d0 *x116*x538 + 0.00224786709554261d0*x117*x595 - 0.0022
     &4786709554261d0*x118*x594 - 0.00142842107214391d0*x191*x252 + 0.00
     &0113462377479451d0*x535) + x591 + x592 - x593 + 0.0501389896966688
     &d0*x261*x535/x121**3 - 1.12219429262413d0*x535/ (x121**5*x260**2))
     &
      x597 = rho_c**(-3.33333333333333d0)
      x598 = x12*x597
      x599 = x290*x9
      x600 = x30*x599
      x601 = 0.121994199821273d0*x29
      x602 = x25*x295*x601
      x603 = x1*x602
      x604 = rho_c**(-5.66666666666667d0)*x1
      x605 = x283*x5
      x606 = x291*(0.0751469112939684d0*x393 - 0.01915695d0*x514 + 0.003
     &94098931294027d0*x598 + 0.00836319657736588d0*x605)
      x607 = rho_c**(-4.33333333333333d0)
      x608 = x1*x135
      x609 = 8.0d0*x43*x63
      x610 = x305**2
      x611 = x271*x610
      x612 = 6.0d0*x303
      x613 = 6.0d0*x304
      x614 = x39**(-1.33333333333333d0)
      x615 = x38*x614
      x616 = x41**(-1.33333333333333d0)
      x617 = x38*x616
      x618 = -x615 - x617
      x619 = x271*x43
      x620 = x342*x619
      x621 = 3.0d0*x45
      x622 = x339*x340*x621
      x623 = x621/x61**2
      x624 = x339*x623
      x625 = 2.16703006958462d0*x538
      x626 = x53*x605
      x627 = 12.0326097596534d0*x626
      x628 = x329*x333
      x629 = 2.476605793811d0*x50*x52
      x630 = rho_c**(-3.16666666666667d0)*x323
      x631 = x629*x630
      x632 = x323*x597
      x633 = x313*x632
      x634 = rho_c**(-3.5d0)*x334
      x635 = x323*x634
      x636 = x314**2
      x637 = x47*x636
      x638 = 9.90642317524398d0*x50
      x639 = rho_c**(-4.16666666666667d0)*x1*x637*x638
      x640 = x1*x312*x59*x607*x636
      x641 = 144.751042466189d0*x636
      x642 = x615 + x617
      x643 = rho_s*(x612 - x613 + x642)
      x644 = rho_c**(-3.16666666666667d0)*x629*x643
      x645 = rho_c**(-3.33333333333333d0)
      x646 = x313*x643*x645
      x647 = x327*x53
      x648 = 2.0d0*x56
      x649 = 33.420578856467d0*x47
      x650 = x312*x649
      x651 = x61/x55**2
      x652 = x647*x651
      x653 = 5.57009647607783d0*x626
      x654 = 2.0d0*x547
      x655 = -x149*x583 - x526*x578 - x549*x654 + 1.94940787421155d0*x55
     &1 - 0.0974703937105774d0*x567 - 20.5190512099364d0*x569 - 0.779763
     &14968462d0*x571 + 20.0d0*x573 - 0.5d0*x574 + x575 + x576 + 5.12976
     &28024841d0*x577 - x584 + x585 - x586 + x587 + x588 - x591 - x592 +
     & x593 + x596
      x656 = 0.803050405822664d0*x154
      x657 = x25*x419
      x658 = rho_c**(-5.0d0)
      x659 = 0.350864882184401d0*x65
      x660 = 24.0d0*x348
      x661 = 48.0d0*x1
      x662 = 2.25322953511185d0*x70
      x663 = x356*x69
      x664 = 0.181712059283214d0*x383
      x665 = 0.0115260005177981d0*x283
      x666 = 0.0201705009061467d0*x645
      x667 = 0.00288150012944953d0*x597
      x668 = x1/x92
      x669 = rho_c**(-4.33333333333333d0)
      x670 = 0.0115260005177981d0*x669
      x671 = 0.817704266774463d0*x21
      x672 = 0.00942485980224818d0*x376
      x673 = 3.63424118566428d0*x89
      x674 = x1*x645
      x675 = x674*x76
      x676 = 3.63424118566428d0*x103
      x677 = x1/x78
      x678 = x362*x503
      x679 = x361*x678
      x680 = 0.0376994392089927d0*x464
      x681 = x316*x76
      x682 = 1.45369647426571d0*x681
      x683 = 0.0582386976490866d0*x393
      x684 = 0.219734258381682d0*x605
      x685 = x684*x94
      x686 = x308*x406
      x687 = 0.116477395298173d0*x20
      x688 = x5*x94
      x689 = x405*x688
      x690 = 0.384534952167943d0*x645
      x691 = 0.0549335645954204d0*x597
      x692 = 0.0582386976490866d0*x20
      x693 = x278*x692
      x694 = 0.219734258381682d0*x669
      x695 = x362*x500
      x696 = x373*x695
      x697 = x363*x389
      x698 = x402*x697
      x699 = x361*x362
      x700 = 0.0141372897033723d0*x21
      x701 = x368*x410
      x702 = x409*x701
      x703 = 1.09027235569928d0*x377
      x704 = x409*x412
      x705 = x684*x82
      x706 = x308*x396
      x707 = x5*x82
      x708 = x395*x707
      x709 = x390*x402
      x710 = 4.36108942279713d0*x384
      x711 = 0.0698864371789039d0*x393
      x712 = 0.329601387572523d0*x5
      x713 = 0.0698864371789039d0*x276
      x714 = 0.329601387572523d0*x288
      x715 = x713*x80 + x714*x82
      x716 = x84**(-3)
      x717 = x391*x716*x88
      x718 = x713*x98 + x714*x94
      x719 = x99**(-3)
      x720 = x391*x411*x719
      x721 = x104*x664 - x362*x373*x700*x702*x92 - x365*x671*x78*(x395*x
     &666 + x395* x667 - x665 - x670*x677) - x366*x672 - x371*x373*x672 
     &+ x371*x671 *(x405*x666 + x405*x667 + x665 + x668*x670) - x381*x68
     &2 + x387* x682 + x39*x709*x710 - x390*x403*x703 - x392*x78*(-x677*
     &x693*x80 - x677*x694*x707 - x683*x80 + x687*x706 + x690*x708 + x69
     &1*x708 - x705) - x403*x717*(-x396*x711 - x398*x712 + x715) - x409*
     &x720*x92 *(x406*x711 + x407*x712 + x718) - x41*x704*x710 + x414*(x
     &668*x688 *x694 + x668*x693*x98 + x683*x98 + x685 + x686*x687 + x68
     &9*x690 + x689*x691) + x664*x90 - x673*x675 - x675*x676 + x679*x680
     & - x680* x696 - x698*x699*x700*x78 - x703*x704*x92
      x722 = 12.0d0*x109
      x723 = x435*x9
      x724 = x157*x723
      x725 = x436*(-x18*x294 + 0.0259624262672375d0*x24*x393 - 0.1519135
     &02485936d0* x24*x605 + 0.455740507457808d0*x288*x434 + 0.037573455
     &6469842d0* x393 + 0.00197049465647014d0*x598 + 0.00418159828868294
     &d0*x605)
      x726 = 12.0d0*x166
      x727 = 1.0525946465532d0*x426
      x728 = 3.61230841029072d0*x174
      x729 = x176*x4*x728
      x730 = x1*x729
      x731 = 5.0d0*x38
      x732 = x40*x731
      x733 = x42*x731
      x734 = 12.0d0*x183
      x735 = 0.0833333333333333d0*rho_c
      x736 = -x216 + x217
      x737 = x146*x736
      x738 = x211*x470
      x739 = 2.56488140124205d0*x216 - 2.56488140124205d0*x217 - x561 - 
     &x563
      x740 = x565*x739
      x741 = x224*x736
      x742 = x224*x550
      x743 = x265*x736
      x744 = x484 - x487
      x745 = x271*x482
      x746 = x344*x736
      x747 = 2.56488140124205d0*x137
      x748 = 4.0d0*x490
      x749 = -x146*x746 + x491 + x572*x748 + x741*x747
      x750 = x408*x412
      x751 = x390*x401
      x752 = x103*x510 - x507*x89 - x511*x751 + x512*x750
      x753 = x505 + x752
      x754 = -x461 + x462
      x755 = x518*x754*x76 + x520 + x752
      x756 = rho_s*x278
      x757 = x271*x305
      x758 = 3.0d0*x304
      x759 = 3.0d0*x303
      x760 = x305*x56
      x761 = x486*x760
      x762 = 1.2383028969055d0*x305*x538
      x763 = x618 + x758 - x759
      x764 = 2.476605793811d0*x321
      x765 = x763*x764
      x766 = x305*x52
      x767 = x630*x638*x766
      x768 = x53*x535
      x769 = x763*x768
      x770 = 18.0938803082737d0*x628
      x771 = 14.859634762866d0*x191 + 33.420578856467d0*x310
      x772 = 1.0d0*x213
      x773 = x134*x247*x34*x490 + x219*x772 + x228*x229 - 10.25952560496
     &82d0*x229* x736 + x229*x739 - x246*x746 - 0.0974703937105774d0*x37
     &*x740 - 1.55952629936924d0*x471*x472 - 16.0d0*x471*x572 + x578*x74
     &8 - x737*x772 + x743*x747
      x774 = x31*x362
      x775 = 0.00471242990112409d0*x774
      x776 = 0.0188497196044964d0*rho_c**(-2.66666666666667d0)*rho_s
      x777 = 0.726848237132856d0*x383
      x778 = 0.0188497196044964d0*x376
      x779 = 0.00864450038834858d0*x315
      x780 = 0.00288150012944953d0*x283
      x781 = 0.0115260005177981d0*x645
      x782 = x408*x701
      x783 = x401*x697
      x784 = 0.0582386976490866d0*x393
      x785 = x784*x80
      x786 = 0.164800693786261d0*x535
      x787 = 0.0549335645954204d0*x605
      x788 = 0.219734258381682d0*x645
      x789 = x784*x98
      x790 = -x359*x39*x699*x783 + x362*x374*x41*x782 + x365*x776 + x370
     &*x776 - x378* x39*x751 + x378*x41*x750 + x381*x777 + x385*x750 + x
     &385*x751 - x387*x777 - x39*x392*(x692*x706 + x708*x788 - x785 - x7
     &86*x82 - x787*x82) - x39*x502*x698 - x402*x511*x715*x716*x88 + x40
     &9*x411* x512*x718*x719 - x41*x413*(x686*x692 + x689*x788 + x786*x9
     &4 + x787*x94 + x789) + x41*x502*x702 + x500*x775 - x503*x775 - x50
     &7* x709 + x510*x704 - x671*x678*(x395*x781 - x779 - x780) - x671* 
     &x695*(x405*x781 + x779 + x780) + x673*x681 + x676*x681 - x679* x77
     &8 + x696*x778
      x791 = 1.16964472452693d0*x472
      x792 = x559 + x562
      x793 = 0.0833333333333333d0*x565*x792
      x794 = x133*x792
      x795 = 20.5190512099364d0*x525
      x796 = 12.0d0*x475
      x797 = x614 + x616
      x798 = -0.854960467080683d0*x137*x794 + x2*x791 - x2*x796 - x490*x
     &654 + x490* x795 - x793
      x799 = x764*x797
      x800 = x768*x797
      x801 = 9.90642317524398d0*x322
      x802 = mu*x535
      x803 = mu*x47*x486
      x804 = 0.0282745794067445d0*x774
      x805 = 0.0141372897033723d0*x501
      x806 = 3.63424118566428d0*x383
      x807 = 4.36108942279713d0*x377
      x808 = x103*x806 + x364*x804 + x369*x804 - x392*(x705 + x785) + x4
     &01*x715*x717 + x408*x718*x720 - x413*(x685 + x789) + x750*x807 + x
     &751*x807 + x782*x805 + x783*x805 + x806*x89
      E = -rho_c*(x139 + x150 + x189 + x190)
      d1E(1) = -0.25d0*rho_c*(x148*x247 + x187*x469 - 0.303232673413971d
     &0*x196 - 0.136878972689611d0*x205 - 4.0d0*x210 + 1.55952629936924d
     &0*x215 + x222 - x227*x228 + 16.0d0*x230 - 4.0d0*x266 + 1.356310494
     &81572d0* x274*x277) + x152 - x189
      d1E(2) = rho_c*x521
      d2E(1) = -0.194940787421155d0*x148*x527 - 0.678155247407862d0*x185
     &*x277 - 0.5d0* x187*(x110*(-x26*x421 + x3*x423 - x418*(x350 + x354
     & + x358 + x367 + x375 + x379 + x380 - x386 + x388 + x404 + x415 + 
     &x74) + x422* x532 + x428) - x152*x154*x309 + x156*(-x149*x527 + x2
     &12*x476 + x346 - x474*x524 - x528 - x530 + x531) - x167*(x430 - x4
     &32 + x439 + x441 - x442) - x184*(x445 - x446 + x448 + x450 + x451 
     &- x467 - x468) + x282 - x297*x533 + x299*x533 + x323*x534 + x33*x5
     &32 - x342*x480) + 0.151616336706985d0*x196 - x206*( -0.00193804554
     &230887d0*x208 + x523) - 0.77976314968462d0*x215 - 0.5d0*x222 + x47
     &5*x526 - 0.0684394863448054d0*x522*x6 - 5.1297628024841d0*x524*x52
     &5 + 2.0d0*x531 - x735*( 1.16964472452693d0*x148*x583 + 6.898418093
     &80227d0*x183*x274*x598/ x186**6 + x187*(-rho_c**(-6.0d0)*x603 - rh
     &o_c**(-6.0d0)*x603 + x153*x655*x656 - 4.28293549772087d0*x154*x308
     &*x347 + 1.18970430492246d0*x270*x308 + 7.85204841248826d0*x278*x30
     &7 - 1.07073387443022d0*x279*x600 + 3.39065726902902d0*x280*x604 - 
     &0.892278228691848d0*x29*x292*x298 + x292*x296*x601 - 0.28465313291
     &6304d0*x293*x295*x426 + 0.000542507213303895d0*x293* x607 + 0.2676
     &83468607554d0*x30*x31*x606 + 3.0d0*x300*x341*(rho_c **(-4.5d0)*x1*
     &x45*x58*x641/x43**4 + x338*x647*x648 + x62*(-x625 + x631 + 11.1401
     &929521557d0*x633 - x639 - 33.420578856467d0*x640 - x644 - 5.570096
     &47607783d0*x646 - x653) + x625 + x627 + 20.3556153468079d0*x628 - 
     &x631 - 24.0652195193068d0*x633 + 18.0938803082737d0*x634*x643 - 54
     &.281640924821d0*x635 + x639 + 72.1956585579204d0*x640 + x644 + 12.
     &0326097596534d0*x646 + x652*( 3.71490869071649d0*x191 + 8.35514471
     &411675d0*x310 - x317*x650 - 14.859634762866d0*x324)) + x302*x566*(
     &-x612 + x613 + x618) + 12.0d0*x306*x620 + x56*x619*x622*(-x311 + x
     &319 + x337) + x608* x609*x611 + x619*x624*(x326 + x328 + x331 - x3
     &32 - x336) + x722*( -0.159903339455843d0*rho_c**(-5.33333333333333
     &d0)*x1*x657 + 0.00910930363347549d0*rho_c**(-3.66666666666667d0)*x
     &269 + 1.75432441092201d0*x1*x658*x66 - 0.00495710127051027d0*x106*
     &x308 - x162*x292*x659 + x267*x497 - 0.133252782879869d0*x268*x419*
     &x597 + 0.0799516697279215d0*x292*x420 - 0.0037178259528827d0*x298*
     &x417 + x418*(-4.50645907022371d0*x1*x597*x663*x70 - 0.751076511703
     &951d0*x10*x355*x356 - x162*x348*x661*x71 + x294* x351*x661 - x3*x3
     &48*x357*x662 - 0.154912426780983d0*x31*x355 - x349*(0.917202013581
     &841d0*x393 - 1.99883265211535d0*x605) - x353* x489 - x660*x72 + x7
     &21) + x422*x606 - 0.701729764368802d0*x427* x599*x65) + x726*(0.21
     &1916079314314d0*x159*x604 - 0.074356519057654d0*x160*x308 - 2.1051
     &8929310641d0*x163*x294*x347 + x165*x655 - 0.133841734303777d0*x279
     &*x724 + 3.15778393965961d0* x426*x440 + 0.111534778586481d0*x431*x
     &437 - x438*x725) + x734*( 1.3157433081915d0*x158*x444*x658 + x182*
     &x655 + 4.82889266566159d0 *x272*x598 - 4.13905085628136d0*x284*x34
     &7 + 0.00619637658813783d0 *x294*x449 - x3*x515 - 0.004647282441103
     &38d0*x431*x465 + 0.526297323276602d0*x433*x437 - x444*x723*x727 - 
     &x447*x725 + x466 *(-rho_c**(-3.66666666666667d0)*x730 - 0.299209d0
     &*x162*x459 - 0.181712059283214d0*x172*x382 - x175*x3*x455*(0.62002
     &8561181324d0 *x393 - 0.156676175559749d0*x605) + 2.42282745710952d
     &0*x178*x382 - x3*x452*x460*x728 - x308*x730 + 2.90739294853142d0*x
     &316*x463 - 14.5369647426571d0*x376*x454 - 43.6108942279714d0*x453*
     &x458 + 36.3424118566428d0*x456*x674 + 0.726848237132856d0*x681*(-6
     &.0d0* x461 + 6.0d0*x462 + x732 + x733) + x721))) - 0.0036946291400
     &831d0 *x191*x545 + 0.410636918068833d0*x201*x205 + 96.0d0*x212*x57
     &8 - 6.78155247407862d0*x274*x275*x393*x45 + 8.13786296889434d0*x27
     &7* x469 + 21.8783328460551d0*x536 - 0.578073242843334d0*x537 + 1.0
     &613143569489d0*x539 - 0.136878972689611d0*x541 + 0.011628273253853
     &2d0*x546 + 24.0d0*x547*x549 - 23.3928944905386d0 *x551 - 0.4106369
     &18068833d0*x558*x6 + 1.16964472452693d0*x567 + 246.228614519237d0*
     &x569 + 9.35715779621544d0*x571 - 240.0d0*x573 + 6.0d0*x574 - 12.0d
     &0*x575 - 12.0d0*x576 - 61.5571536298092d0* x577 - 12.0d0*x596)
      d2E(2) = x521 - x735*(1.16964472452693d0*x0*x740 + 3.0d0*x147*x246
     &*x37*x736 + x187*(rho_s*x602*x658 + x0*x302*(x642 - x758 + x759) -
     & 2.14146774886044d0*x154*x278*x749 + 6.0d0*x220*x341*x611*x744 + x
     &271*x622*x761 - 1.96301210312207d0*x280*x756 + x308*x656*x773 + x3
     &23*x489*x609*x757 + 0.535366937215109d0*x477*x600 + 6.0d0*x479* x6
     &20 - x484*x624*x757 + x722*(rho_s*x294*x599*x659 + 0.0799516697279
     &215d0*rho_s*x607*x657 - rho_s*x66*x727 - x27*x424* x753 + x418*(rh
     &o_s*x15*x660*x71 - 48.0d0*rho_s*x162*x351 + rho_s* x283*x662*x663 
     &+ x220*x352 + x790)) + x726*(-0.122688256445129d0* x159*x756 - x16
     &3*x727*x749 + 0.0669208671518886d0*x477*x724 + x496*x773) + x734*(
     &-0.789445984914903d0*rho_s*x158*x17*x426 + 0.526297323276602d0*rho
     &_s*x514*x723 - 0.00232364122055169d0*x431* x755 + x466*(rho_s*x153
     &*x729 - 29.0739294853142d0*x316*x456 + 1.45369647426571d0*x383*x75
     &4 + 0.726848237132856d0*x383*(3.0d0* x461 - 3.0d0*x462 - x732 - x7
     &33) + 21.8054471139857d0*x453*x464 + x790) + x517*x773 - 2.0695254
     &2814068d0*x598*x749) + 3.0d0*x745*( -x305*x627 - 27.1408204624105d
     &0*x305*x628 + 72.1956585579204d0* x305*x633 + 144.751042466189d0*x
     &305*x635 - x305*x652*x771 - x338* x53*x761 + x484*x647*x760 + x62*
     &(-x305*x632*x650 + x305*x653 + x762 - x765 - x767 - 5.570096476077
     &83d0*x769) - x762 + x763*x770 + x765 + x767 + 12.0326097596534d0*x
     &769)) - 30.7785768149046d0* x212*x743 - 123.114307259618d0*x218*x7
     &42 + 4.06893148444717d0* x277*(x110*(x418*x753 + x498) + x167*(x49
     &5 + x496*x749) + x184*( x466*x755 + x516 + x517*x749) + x302*x479 
     &+ x305*x744*x745 + x478 + x493*x749) - 4.67857889810772d0*x34*x471
     &*x570 + 12.0d0*x34*x550 *x737 - 48.0d0*x471*x578 + 18.714315592430
     &9d0*x472*x738 - 12.0d0* x547*x550 + 123.114307259618d0*x550*x741 +
     & 192.0d0*x572*x738 - 12.0d0*x739*x742)
      d2E(3) = -rho_c*(x0*x793 + 0.0833333333333333d0*x187*(-x0*x36*x609
     &*x636 + x278* x656*x798 + 0.535366937215109d0*x280*x308 - x314*x34
     &0*x486*x488* x52*x56*x621 + 12.0d0*x341*x479*x488 - 3.0d0*x36*x482
     &*(x305*x485* x648*x803 + x335*x641 + x610*x651*x771*x803 - x62*(x6
     &10*x649*x802 + x610*x801 + x799 + 5.57009647607783d0*x800) + x636*
     &x801 + 72.1956585579204d0*x637*x802 + x770*x797 + x799 + 12.032609
     &7596534d0*x800) - x484*x488*x623*x766 + x534*x797 + x722 *(-x418*(
     &-x499 + x808) + x497) + x726*(x164*x426*x798 + x494) + x734*(-x466
     &*(-x43*x806 - x519 + x808) + x515 + 1.55214407110551d0 *x598*x798)
     &) + 0.854960467080682d0*x550*x794 - x608*x791 + x608* x796 + x654*
     &x738 - x738*x795)
      end subroutine


C*****************************************************************************
      pure subroutine D2ESRC_VWN5_ERF_singletref_triplet(rho_c, mu, E, d
     &1E, d2E)
C*****************************************************************************
C   Implemented by E.R. Kjellgren.
C
C   Subroutine generated using Sympy 1.3
C   Generated: March 25, 2019
C*****************************************************************************
      implicit none
      real*8, intent(in) :: rho_c, mu
      real*8, intent(out) :: E, d1E(9), d2E(45)
      real*8 :: x0, x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12, x
     &13, x14, x15, x16, x17, x18, x19, x20, x21, x22, x23, x24, x25, x2
     &6, x27, x28, x29, x30, x31, x32, x33, x34, x35, x36, x37, x38, x39
     &, x40, x41, x42, x43, x44, x45, x46, x47, x48, x49, x50, x51, x52,
     & x53, x54, x55, x56, x57, x58, x59, x60, x61, x62, x63, x64, x65, 
     &x66, x67, x68, x69, x70, x71, x72, x73, x74, x75, x76, x77, x78, x
     &79, x80, x81, x82, x83, x84, x85, x86, x87, x88, x89, x90, x91, x9
     &2, x93, x94, x95, x96, x97, x98, x99, x100, x101, x102, x103, x104
     &, x105, x106, x107, x108, x109, x110, x111, x112, x113, x114, x115
     &, x116, x117, x118, x119, x120, x121, x122, x123, x124, x125, x126
     &, x127, x128, x129, x130, x131, x132, x133, x134, x135, x136, x137
     &, x138, x139, x140, x141, x142, x143, x144, x145, x146, x147, x148
     &, x149, x150, x151, x152, x153, x154, x155, x156, x157, x158, x159
     &, x160, x161, x162, x163, x164, x165, x166, x167, x168, x169, x170
     &, x171, x172, x173, x174, x175, x176, x177, x178, x179, x180, x181
     &, x182, x183, x184, x185, x186, x187, x188, x189, x190, x191, x192
     &, x193, x194, x195, x196, x197, x198, x199, x200, x201, x202, x203
     &, x204, x205, x206, x207, x208, x209, x210, x211, x212, x213, x214
     &, x215, x216, x217, x218, x219, x220, x221, x222, x223, x224, x225
     &, x226, x227, x228, x229, x230, x231, x232, x233, x234, x235, x236
     &, x237, x238, x239, x240, x241, x242, x243, x244, x245, x246, x247
     &, x248, x249, x250, x251, x252, x253, x254, x255, x256, x257, x258
     &, x259, x260, x261, x262, x263, x264, x265, x266, x267, x268, x269
     &, x270, x271, x272, x273, x274, x275, x276, x277, x278, x279
      E = 0.0d0
      d1E(:) = 0.0d0
      d2E(:) = 0.0d0
      x0 = mu**2
      x1 = 2.14502939711103d0
      x2 = 1d0/x1
      x3 = rho_c**0.666666666666667d0
      x4 = 1d0/x3
      x5 = x2*x4
      x6 = 0.508616435555896d0*x0*x5 + 1.0d0
      x7 = x6**(-4)
      x8 = mu**4
      x9 = 1.46459188756152d0
      x10 = 1d0/x9
      x11 = rho_c**0.333333333333333d0
      x12 = 1d0/x11
      x13 = x10*x12
      x14 = exp(-0.6836107611867116d0*x13)
      x15 = rho_c**(-1.33333333333333d0)
      x16 = 4.60115111447049d0
      x17 = 1d0/x16
      x18 = x15*x17
      x19 = rho_c**1.0d0
      x20 = 1d0/x19
      x21 = 3.14159265358979d0
      x22 = 1d0/x21
      x23 = 0.009578475d0*x22
      x24 = 0.0188171922990732d0*x13 + 0.00126674656487366d0*x18 - x20*x
     &23 + 0.0676322201645715d0*x5 + 1.0d0
      x25 = x14*x24
      x26 = x25 - 1.0d0
      x27 = 0.13157433081915d0*x26
      x28 = x20*x22
      x29 = x27*x28
      x30 = 2.28942848510666d0
      x31 = x16*x30
      x32 = 1.63540853354893d0*x31
      x33 = exp(-0.49698248213959023d0*x13)
      x34 = x33*(-0.352521395009435d0*x13 + 0.558025705063192d0*x5)
      x35 = x1*x34
      x36 = 10.9027235569928d0*x35
      x37 = 0.494402081358784d0*x13 + 0.0524148278841779d0*x5 + 1.0d0
      x38 = 1d0/x37
      x39 = -0.0259335011650457d0*x13 + 1.0d0
      x40 = x38*x39
      x41 = 1.63540853354893d0*x31
      x42 = x40*x41
      x43 = 0.148394183600699d0
      x44 = 0.00139418473233101d0*rho_c**(-1.0d0)*x43
      x45 = x44*(-x32 + x36 + x42)
      x46 = 0.826307487110758d0
      x47 = rho_c**(-0.166666666666667d0)*x46
      x48 = 1.90636858599387d0*x47
      x49 = x48 + 3.72744d0
      x50 = 0.90856029641607d0*x13
      x51 = 3.5529372610885d0*x47 + x50 + 12.9352d0
      x52 = 1d0/x51
      x53 = 0.953184292996937d0*x47
      x54 = x53 + 0.10498d0
      x55 = x54**2
      x56 = x52*x55
      x57 = 0.000969022771154437d0*log(x56) + 0.0310907d0*log(x50*x52) +
     & 0.038783294878113d0*atan(6.1519908197590798d0/x49)
      x58 = 1.55214407110551d0*x18
      x59 = x57*x58
      x60 = 0.101321183642338d0
      x61 = 2.0d0*x60
      x62 = 0.693147180559945d0
      x63 = -x62 + 1.0d0
      x64 = x0*x13
      x65 = 5.57236303607474d0*mu*x47 + 1.0d0
      x66 = 3.13317926779378d0*x64 + x65
      x67 = 1d0/x66
      x68 = 0.564189583547756d0
      x69 = mu**3
      x70 = 3.39260255780131d0*rho_c**(-0.5d0)*x68*x69 + 6.7683429898050
     &4d0*x64 + x65
      x71 = x67*x70
      x72 = log(x71)
      x73 = x63*x72
      x74 = mu**6
      x75 = rho_c**(-1.66666666666667d0)
      x76 = x43*x75
      x77 = 0.101321183642338d0
      x78 = x57*x77
      x79 = rho_c**(-2.0d0)
      x80 = 0.526297323276602d0*x79
      x81 = 0.179587122125167d0
      x82 = x25*x81
      x83 = 0.175432441092201d0*x20
      x84 = 0.0837224526465878d0
      x85 = -4.49737346725955d0*x13 + 0.825481812223657d0*x5
      x86 = exp(-0.28165369188898165d0*x13)
      x87 = x21*x86
      x88 = 12.0d0*x87
      x89 = x85*x88
      x90 = x84*(x19*x89 + x3*x42)
      x91 = 0.00111534778586481d0*x75
      x92 = 1.41421356237310d0
      x93 = x69*x92
      x94 = mu**8
      x95 = 0.0472353356922751d0
      x96 = x57*x95
      x97 = rho_c**(-2.66666666666667d0)
      x98 = 0.0669208671518886d0*x97
      x99 = 0.0223069557172962d0*x75
      x100 = mu**5
      x101 = x25*x84
      x102 = x101*x92
      x103 = x100*x102
      x104 = -x103*x99 + x74*(-0.0167302167879722d0*x26*x76 + x78*x80) -
     & x93*(x82*x83 + x90*x91) + x94*x96*x98
      x105 = x104 - x61*x73
      x106 = x57 - x7*(x105 - x8*(x29 + x45 - x59))
      x107 = rho_c**(-1.16666666666667d0)*x46
      x108 = x49**(-2)
      x109 = 37.8469910464d0*x108 + 1.0d0
      x110 = 1d0/x109
      x111 = x108*x110
      x112 = x107*x111
      x113 = 0.0758081683534927d0*x112
      x114 = 0.30285343213869d0*x15
      x115 = x10*x114 + 0.592156210181417d0*x107
      x116 = x115*x52
      x117 = 0.90856029641607d0*x12
      x118 = x116*x117
      x119 = x114 - x118
      x120 = 0.0342197431724027d0*x11
      x121 = x119*x120
      x122 = 1d0/x54
      x123 = 0.000307885761673592d0*x107
      x124 = x116*x54
      x125 = 0.000969022771154437d0*x124
      x126 = x122*(x123 - x125)
      x127 = -x29 + x59
      x128 = x6**(-5)
      x129 = x2*x75
      x130 = x128*x129
      x131 = x0*x130*(x105 + x8*(x127 - x45))
      x132 = x22*x79
      x133 = x132*x27
      x134 = x32 - x36 - x42
      x135 = x134*x43
      x136 = 0.00232364122055169d0*x135*x79
      x137 = x10*x15
      x138 = rho_c**(-2.33333333333333d0)
      x139 = x138*x17
      x140 = -0.0225440733881905d0*x129 + 0.0047892375d0*x132 - 0.003136
     &1987165122d0* x137 - 0.000844497709915773d0*x139
      x141 = x14*(0.113935126864452d0*x137*x24 + x140)
      x142 = 0.263148661638301d0*x28
      x143 = x141*x142
      x144 = 2.06952542814068d0*x139*x57
      x145 = -x114 + x118
      x146 = x120*x145
      x147 = x113 + x122*(-x123 + x125) + x146
      x148 = x147*x58
      x149 = rho_c**(-0.333333333333333d0)
      x150 = x149*x31
      x151 = x4*x9
      x152 = x33*(-0.372017136708795d0*x129 + 0.117507131669812d0*x137)
      x153 = x1*x152
      x154 = 10.9027235569928d0*x3
      x155 = x21*x30*x38
      x156 = 0.0141372897033723d0*x155*x4
      x157 = 1.09027235569928d0*x150*x40
      x158 = x37**(-2)
      x159 = 0.034943218589452d0*x129 + 0.164800693786261d0*x137
      x160 = x158*x159
      x161 = x160*x39
      x162 = x3*x41
      x163 = x161*x162
      x164 = -7.26848237132856d0*x149*x35 + 1.09027235569928d0*x150 - 1.
     &80615420514536d0*x151*x34 - x153*x154 - x156 - x157 - x163
      x165 = 0.00139418473233101d0*x76
      x166 = x164*x165
      x167 = x43*x97
      x168 = 0.0278836946466203d0*x167*x26
      x169 = 0.0334604335759443d0*x76
      x170 = x141*x169
      x171 = rho_c**(-3.0d0)
      x172 = 1.0525946465532d0*x171*x78
      x173 = x77*x80
      x174 = x147*x173
      x175 = 0.00168899541983155d0*x139
      x176 = x23*x79
      x177 = 0.045088146776381d0*x129
      x178 = 0.00627239743302441d0*x137
      x179 = x14*(-x175 + x176 - x177 - x178)
      x180 = x179*x81
      x181 = x180*x83
      x182 = 0.122619224952946d0
      x183 = x182*x25
      x184 = 0.0399758348639607d0*x138*x183
      x185 = 0.175432441092201d0*x79*x82
      x186 = 0.00185891297644135d0*x90*x97
      x187 = x85*x86
      x188 = 2.14502939711103d0
      x189 = rho_c**(-0.333333333333333d0)*x188
      x190 = -0.550321208149104d0*x129 + 1.49912448908652d0*x137
      x191 = x19*x88
      x192 = x84*(x156 + x157 + x163 + 1.12661476755593d0*x187*x189 + x1
     &90*x191 + x89 )
      x193 = x192*x91
      x194 = x0*x92
      x195 = 1d0/x70
      x196 = 0.928727172679123d0*x107
      x197 = mu*x137
      x198 = 2.25611432993501d0*x197
      x199 = x0*x68
      x200 = 1.69630127890066d0*rho_c**(-1.5d0)*x199
      x201 = x196 + 1.04439308926459d0*x197
      x202 = -x196 - x198 - x200 + x201*x71
      x203 = x195*x202
      x204 = mu**7
      x205 = x204*x95
      x206 = x205*x98
      x207 = x8*x92
      x208 = x207*x84
      x209 = x208*x99
      x210 = rho_c**(-3.66666666666667d0)
      x211 = x204*x96
      x212 = x102*x8
      x213 = 0.0571643564037363d0
      x214 = x207*x25
      x215 = x213*x214
      x216 = -0.00508309165921971d0*x171*x215 - 0.17845564573837d0*x210*
     &x211 + 0.037178259528827d0*x212*x97
      x217 = mu*x7
      x218 = x217*(-x100*(-x168 + x170 + x172 - x174) + x14*x209*(x175 -
     & x176 + x177 + x178) - x194*(x181 + x184 - x185 - x186 + x193) - x
     &203*x61*x63 + x206*(x113 - x121 - x126) + x216 - x69*(-x133 + x136
     & + x143 + x144 - x148 - x166))
      x219 = 0.000615771523347183d0*x107
      x220 = rho_c**(-2.33333333333333d0)*x10
      x221 = 1.82319440383792d0*x220/(x109**2*x49**5)
      x222 = 0.0481727702369445d0*x110*x220/x49**3
      x223 = rho_c**(-2.16666666666667d0)*x46
      x224 = 0.0884428630790749d0*x111*x223
      x225 = 0.0114065810574676d0*rho_c**(-0.666666666666667d0)*x145
      x226 = 1d0/x55
      x227 = -0.317728097665645d0*x107 + x124
      x228 = x123*x226*x227
      x229 = x116*x146
      x230 = 0.000969022771154437d0*x116*x122*x227
      x231 = 0.40380457618492d0*x138
      x232 = x10*x231 + 0.69084891187832d0*x223
      x233 = x115*(1.18431242036283d0*x107 + 0.60570686427738d0*x137)/x5
     &1**2
      x234 = x120*(-0.60570686427738d0*x116*x15 - x117*x232*x52 + x117*x
     &233 + x231)
      x235 = x226*(-x124*x219 + 4.89119786774443d-5*x220 + 0.00035920005
     &5285857d0* x223*x54 - 0.000969022771154437d0*x232*x56 + 0.00096902
     &2771154437d0*x233*x55)
      x236 = x61*(x62 - 1.0d0)
      x237 = x104 + x236*x72 + x8*(x127 + x134*x44)
      x238 = rho_c**(-3.33333333333333d0)
      x239 = x17*x238
      x240 = x2*x97
      x241 = x203*x236
      x242 = x171*x22
      x243 = x10*x138
      x244 = x14*(0.00394098931294027d0*x239 + 0.0751469112939684d0*x240
     & - 0.01915695d0*x242 + 0.00836319657736588d0*x243)
      x245 = rho_c**(-4.0d0)
      x246 = rho_c**(-4.66666666666667d0)
      x247 = x196 + x198 + x200
      x248 = mu*x201
      x249 = x248*x67
      x250 = 1.08351503479231d0*x223
      x251 = mu*x243
      x252 = rho_c**(-2.5d0)*x199
      x253 = -x221 + x222 - x224 + x225 + x228 - x229 - x230 + x234 + x2
     &35
      x254 = 0.350864882184401d0*x171*x82
      x255 = 24.0d0*x87
      x256 = rho_c**(-1.33333333333333d0)
      x257 = x256*x31
      x258 = 0.0282745794067445d0*x160*x21*x30*x4
      x259 = 0.219734258381682d0*x243
      x260 = x162*x39
      x261 = x158*x260
      x262 = x150*x161
      x263 = x159*x260*(0.0698864371789039d0*x129 + 0.329601387572523d0*
     &x137)/x37**3
      x264 = 0.363424118566428d0*x257*x40 - x258 + x261*(0.0582386976490
     &866d0*x240 + x259) - 2.18054471139857d0*x262 - x263
      x265 = x84*x91
      x266 = x210*x26*x43
      x267 = x14*(0.455740507457808d0*x137*x140 - x171*x23 + 0.001970494
     &65647014d0* x239 + 0.0259624262672375d0*x24*x240 - 0.1519135024859
     &36d0*x24* x243 + 0.0375734556469842d0*x240 + 0.00418159828868294d0
     &*x243)
      x268 = 0.263148661638301d0*x242*x26
      x269 = x256*x35
      x270 = atan(7.1231089178181177d0/(x48 + 1.13107d0))
      x271 = 1d0/(1.07811815828004d0*x47 + x50 + 13.0045d0)
      x272 = log(x271*x50)
      x273 = log(x271*(x53 + 0.0047584000000000003d0)**2)
      x274 = x60/rho_c**2
      x275 = 0.317708004743941d0*x270 + x272 + 0.000414033794282063d0*x2
     &73
      x276 = 1.2383028969055d0*rho_c**(-2.16666666666667d0)*x46
      x277 = mu*x220
      x278 = 7.26848237132856d0*x257
      x279 = 0.0565491588134891d0*x155*x75 + x258 - x261*(0.058238697649
     &0866d0*x240 + x259) + 8.72217884559427d0*x262 + x263 + x278*x40
      E = rho_c*x106
      d1E(1) = -rho_c*(-x113 + x121 + x126 + 1.35631049481572d0*x131 + x
     &218) + x106
      d2E(1) = -rho_c*(-2.26051749135954d0*x0*x128*x237*x240 + 2.7126209
     &8963145d0*x130* x69*(x100*(x168 - x170 - x172 + x174) + x147*x206 
     &- x179*x209 + x194*(-x181 - x184 + x185 + x186 - x193) + x216 + x2
     &41 + x69*( x133 - x136 - x143 - x144 + x148 + x166)) + x217*(mu*x2
     &02*x236* x247/x70**2 - 4.52089344419912d-5*rho_c**(-4.333333333333
     &33d0)* x214 + x100*(0.111534778586481d0*x141*x167 - 2.105189293106
     &41d0* x147*x171*x77 - x169*x267 + x173*x253 + 3.15778393965961d0*x
     &245* x78 - 0.074356519057654d0*x266) - 0.356911291476739d0*x147*x2
     &05* x210 - 0.0101661833184394d0*x171*x179*x207*x213 + 0.0743565190
     &57654d0*x179*x208*x97 + x194*(-0.00910930363347549d0* rho_c**(-3.6
     &6666666666667d0)*x101 - 0.0799516697279215d0*x138* x179*x182 + 0.3
     &50864882184401d0*x180*x79 + 0.133252782879869d0* x183*x238 + 0.003
     &7178259528827d0*x192*x97 - 0.00495710127051027d0 *x210*x90 - x244*
     &x81*x83 - x254 + x265*(-0.751076511703951d0*x15* x187*x188 - 0.154
     &912426780983d0*x187*x75 - 2.25322953511185d0* x189*x190*x86 - x190
     &*x255 - x191*(0.917202013581841d0*x240 - 1.99883265211535d0*x243) 
     &+ x264)) + x195*x236*(-2.0d0*x247*x249 + x248*x70*(1.8574543453582
     &5d0*x107 + 2.08878617852919d0*x197)/x66 **2 + x250 + 3.00815243991
     &335d0*x251 + 2.54445191835098d0*x252 - x71*(x250 + 1.3925241190194
     &6d0*x251)) + x206*x253 - x209*x244 - 0.0991420254102054d0*x210*x21
     &2 + 0.654337367707355d0*x211*x246 + 0.023721094409692d0*x215*x245 
     &- x241*x249 + x69*( 0.526297323276602d0*x132*x141 + 0.006196376588
     &13783d0*x135*x171 - 4.13905085628136d0*x139*x147 - x142*x267 - 0.0
     &0464728244110338d0* x164*x167 + x165*(-x1*x154*x33*(0.620028561181
     &324d0*x240 - 0.156676175559749d0*x243) - 14.5369647426571d0*x149*x
     &153 - 3.61230841029072d0*x151*x152 - 0.363424118566428d0*x257 + x2
     &64 + 2.42282745710952d0*x269 - 0.299209d0*x34*x79) + 4.82889266566
     &159d0*x239*x57 + x253*x58 - x268)) + x221 - x222 + x224 - x225 - x
     &228 + x229 + x230 - x234 - x235 + 2.29947269793409d0*x237*x239*x8/
     &x6**6) - 0.0684394863448054d0*x11 *x119 + 0.151616336706985d0*x112
     & - x122*(-0.00193804554230887d0* x124 + x219) - 2.71262098963145d0
     &*x131 - 2.0d0*x218
      d2E(3) = -rho_c*(x274*(0.0529513341239902d0*x270 + 0.1666666666666
     &67d0*x272 + 6.90056323803438d-5*x273) + 0.333333333333333d0*x7*(-6
     &.0d0*mu* x195*x60*x63*(2.26173503853421d0*x252 + x276 + 3.00815243
     &991335d0 *x277 - x71*(x276 + 1.39252411901946d0*x277)) + 0.1338417
     &34303777d0*x103*x210 - 0.000160139631553455d0*x246*x275* x94 + 4.0
     &d0*x274*x73 - 3.0d0*x74*(0.000900493163574243d0*x245* x275 - 0.033
     &4604335759443d0*x266) + 3.0d0*x8*(-x165*( -21.8054471139857d0*x269
     & - x278 + x279) - 0.00569658441823587d0* x238*x275 + x268) + 3.0d0
     &*x93*(x254 - x265*(-x20*x255*x85 + x279 ))))
      end subroutine
